Discretizations of a Perturbed Logistic Equation *

The logistic equation has been used to model a population where the intrinsic rate of growth is a linearly decreasing function of the population density. We propose some models arising from discretizations of perturbed logistic equation where external factors such as harvesting 
and migration are included. The different discretizations exhibit diverse dynamic asymptotic behavior when dependence on two control parameters are allowed. We also give an alternative method to provide with new models that produce better approximating solutions.


INTRODUCTION
A general setting for population models is the relation dP(t) dt Births-Deaths :t: Migration where P(t) represents the total population at time t.
Assuming analytical laws depending on the population density for the birth B(P) and death D(P) rates and analytical migration M(P) of the form, B(P) Z akP: D(P) Z/3:Pk k=O k=O M(V) k=O we obtain the differential equation dP dt Z(a -/3 + %)P k=0 aoP bo P2 -{-co P3 q--mo n t-O(p4).(1) This general model includes the well known models of Malthus (with bo, Co, and mo equal to zero) and the logistic equation of Verlhust (Co and mo equal to zero).Here ao and bo are the so-called vital coefficients of the population; aoP is the term of pure exponential growth, ao is known as the per capita growth rate and -bo P2 the term of competi- tion.This last term is justified because individual members of the population compete for limited resources, e.g.land and food.We will consider only the cases where the right-hand side of (1) is a cubic Partially supported by CONACYT.
Corresponding author.126 M.A. MORELES AND F. SOLIS polynomial and the coefficients Co and m0 are small enough, so that the equation might be viewed as a structural perturbation of the logistic equation.
The perturbation terms can model some biological important quantities as harvesting, immigration, emigration, etc.We are considering only the cases where the perturbation terms are simple enough so that the harvest (emigration, immigration, etc) rate is a constant value independent of the population size or it can vary proportional to the third power of the population size.Linear and quadratic powers are already included in the model.
A standard technique to use the continuous sys- tem (1) to study the asymptotic behavior of some species that live in isolated generations is to reduce it to a discrete map.For example, for the logistic equation the following two models have been con- sidered as its difference equation analogue (see [7]), Pn+l (1 + a)Pn(1 Pn) and Pn+, Pn exp(a-(1 q-a)Pn). (2) It is known that there is an optimal discretiza- tion with zero local truncation error, unfortunately there is no procedure to obtain it for a general sys- tem.For this reason we use, as a benchmark, the perturbed logistic equation for which in some cases the optimal discretization can be computed in an explicit form.We use this optimal discretisation to compare it with some discretizations arising from numerical one-step explicit methods for solving dif- ferential equations, e.g.Runge-Kutta methods.In particular models in (2) can be obtained by using a first order Runge-Kutta method.
A drawback of these two models is that their asymptotic behavior coincides with the continuous system only up to some value of a. Moreover any convex combination of these two schemes will not extend the range of validity of the asymptotic solu- tion for larger values of a.To improve the range of validity we will introduce new schemes provided by nonlinear convex combination of Runge-Kutta methods.Of special interest is to study the biologi- cal implications of the discrete maps obtained when we allow dependence on two control parameters.

DISCRETE MODELS
The point of departure is the perturbed logistic equation dP aoP boP 2 + coP -+-too, to simplify our analysis, let y-bo/(ao + 1)P, thus we are led to a equation of the form dy dt ay-(1 + a)y 2 + cy + rn. (4) Euler's method yields the corresponding difference equation Y,+l (1 + a)yn(1 y,) + cy3, + m (5)   which we shall refer as the perturbed free growth model equation.An alternative discretization is obtained by rewriting Eq. ( 4) in the form: dlny dt a-(1 + a)y + Cy 2 q-my -1 and applying the Euler scheme to the resulting equation to obtain: y+l yn exp a (1 + a)y + cy2 + (7) hereafter referred as the perturbed free growth ex- ponential model.
In essence (5) and ( 7) are obtained by applying a one-step numerical method to solve (4) with step size equal to one.Notice that ( 5) and ( 7) are pertur- bations of the systems given in (2).
In regards to truncation error, Euler's methods is a first order method, meaning that the error growths as O(1).It is natural to do the analogue for better numerical methods, like Runge-Kutta methods of higher order (this has been done in [1] for the logistic equation).To illustrate this fact, we consider a second order Runge-Kutta method known as the improved Euler's method.Applying this scheme to model (4) we obtain the following system: z, +l Az, ( + A) + e(rn, A) where A (a2/2) + a + 1, zn fix,, /3 (1 + 2.5a + 2a 2 + 0.5a3)/A, C(x, A) ((1 + a)3x(1 x/2))/ , h(x, A) + (a/2) (1 + a)2x + (1 + a)2x and e(rn, c) a small parameter depending on m and c.System (8) can be considered as a perturbed discrete logistic equation with an extra term that can be interpreted as a regulatory mechanism of harvest- ing or migration.
Our next step is to study numerically the discrete systems ( 5), ( 7) and (8) in order to compare the asymptotic behavior of these models and the effect of discretization.

NUMERICAL RESULTS
For the perturbed logistic equation the behavior of asymptotic solutions are well known but the dynamics of their discretized counterparts are very difficult to analyze (specially for higher order methods).We will give in the next section some description of the dynamics of the discrete models that we propose above.

Constant Migration
In order to study the behavior of the systems ( 5) and (7) under constant migration, we set the value of c to be equal to zero.For system (5), the case m 0 corresponds to the quadratic logistic equation model that has been studied extensively (see [3]).
This case shows a typical period doubling sequence to chaos.The case m > 0 have a similar behavior but the system becomes chaotic for smaller values of the growth rate.The figures that follow are bifurcation diagrams of the grow rate a against population densities.In these diagrams the first thousand points are discarded to let a stable period (or chaos) becomes established and then we plot the next hundred points for each value of a.The ; !: !! , ;; +i',:,' " ;i:: :ii:::-!:::i!!iii.:.":". . . . . . . . . . . . . .: : : : .: ' . . . . .,:..
numerical simulations were carried out on a Sparc station 5.The bifurcation diagrams move up and left with m.This is shown in Fig. 1, where we present several bifurcation diagrams of the system (5), with values of m taken as 0.08, 0.05, 0.0287, 0.0.For m > 0 the diagrams are truncated at the first bifur- cation point, only the case m-0 is shown in full.
In contrast with the previous model which has a homogeneous behavior, the exponential model pre- sents three different characteristics.For values of m < 0.04 the behavior of the system is nonchaotic, the model starts with a single period limit followed by a double period before extinction; for values of rn (0.0287, 0.4) there is a period doubling reversal feature again without chaos (see [9]).Finally for rn < 0.0287 the behavior of the system becomes chaotic in a small interval of the rate growth.The chaotic behavior for very small values of rn presents a similar behavior as model ( 5), but the chaotic region starts for larger values of the parameter a.
In model (5) chaos appears approximately at 2.57 whereas in model ( 7) it appears well after 2.6.In Fig. 2 we present several bifurcation diagrams of the exponential model, with the values of m taken as 0.08, 0.06, 0.04, 0.03.In this case all diagrams are shown completely.
In Fig. 3 we show how chaos is concentrated in a small region of parameter space, in this case we ---2----,.--" .''" --""" took the value of m--0.0287.As m gets smaller the chaotic regime increases and the feature of period doubling reversal disappears, these facts are shown in Figs. 4 and 5 where we set the value of m equal to zero.

Nonconstant Migration
In order to study the behavior of the systems ( 5) and ( 7) under variable (cubic) migration, we set the value of m equal to zero.Consider now the system (5) with m 0 and c > 0. We obtain a similar behavior as in system (5) with c 0 and m > 0 of a period doubling sequence to chaos.The appear- ance of chaos is controlled by the value of c.For larger values of c larger values of a are required for chaos to appear.This effect is shown in Fig. 5 where there are three diagrams corresponding to values of c equal to 0.69, 0.7 and 0.8.Again we only show the case c--0.69 completely and the other two are truncated at the first bifurcation point.In this figure we notice how the value of c delays the appearance of chaos.Another feature to note is the curvature of the one period branches due to the nonlinear perturbation.Compare this effect with those given previously with c 0 and m > 0.
For the exponential model with variable migra- tion, chaos is always present, in contrast with the case of constant migration.Even though chaos is present in this model, its appearance starts far beyond than its counterpart (free growth model with variable migration).Thus the value of c plays the role of a delay of the appearance of chaos in both models, but it has a stronger effect when considering the model (7).In Fig. 6 we show three diagrams for the same values of c as in Fig. 5.Here we note that the chaotic behavior starts for values of a close to four whereas in its counterpart starts approximately at a-2.5.

Migration with Control
The bifurcation diagrams shown in Fig. 7 corre- spond to the system (8) for different values of e(rn, c).We show the complete diagram when e 0, and truncated diagrams at the first bifurcation point when e > 0. The dynamics of this system for e 0 are similar to those when e > 0.
The system exhibits the same characteristic of becoming chaotic for smaller values of the param- eter a when larger values of e are chosen.Recall that the same effect was present in the case of free growth.In this case, we observe three important new features.
(1) The presence of a tangent bifurcation for all values of m.
(2) There is a parameter region where period doubling reversal is present, a feature that has not been documented for polynomial maps before.
(3) Finally and more importantly is the fact that chaos appears for larger values of a than the cases when Euler's method was used.

ALTERNATIVE METHOD OF DISCRETIZATION
From our previous numerical experiments we real- ize that the discretizations obtained are only valid up to some value of the natural parameter.The same feature is present when we use higher order Runge-Kutta methods as shown by the methods above.
As we remarked before, a basic principle to ob- tain a discrete model of population dynamics from a system of differential equations is to apply a nu- merical method.The discrete map need not come from a single numerical method.It is natural to pose the follow- ing problem.
Given a finite collection of numerical methods find a discrete map, better in the sense it extends the range of validity of asymptotic solutions.
We propose a solution as follows.
Let C be the family of continuous functions from R into (0, 1) and let {..}/N__ be a finite collection of maps obtained by discretizations of one given dif- ferential equation.Choose /3; E C for i= 1,...,n with the property that }-/3;(A)= 1. Associate the discrete map given by (9) The choice of the weights (/3's) may depend on the particular differential equation.For instance if we use the weights as constant functions the range of validity improves at least in one hundred percent.In particular with the two Runge-Kutta methods that we use in the previous section and choosing/31 =/3 (the weight for the improved Euler's method) as the value that optimizes the function E(/3)= At., where At. is the value of the parameter where the discre- tization has its first bifurcation, a 160% improve- ment is achieved.We obtain in this case/3 =0.255 and A.= 5.2 for m =0 and c=0.Compare these values with the ones given in Fig. 7.

CONCLUSIONS
We have consid.eredthe perturbed logistic equa- tion (4) to study population evolution with the aid of some discrete models associated with it.Ecological external factors, such as migration, were pre- sent as perturbation of the logistic equation.Such perturbations terms resulted by considering a higher order polynomial for the birth and death rates as in (1) or by adding to the logistic equation itself the corresponding terms.We clarified that two of the most popular discrete models associated with (2) are nothing but different discretizations of the logistic equation.Thus we propose that if external factors are to be considered, they should be added to the continuous equation and then discretize the resulting equation with appropriate methods.
Compare with [7] where the migration term is added to (2).
The dynamics exhibited by the discrete models were very diverse.Phenomena like nonchaotic period doubling and period doubling sequence to chaos, chaotic and nonchaotic period doubling reversal appeared.When chaos is present, we have shown at least numerically that its appearance depends on the method of discretization.More precisely, for better methods, the interval of stabil- ity is larger and the appearance of chaos is delayed.We noticed that in some cases, chaos in the logistic equation can be controlled by appropriate per- turbations; in ecological terms, populations with chaotic behavior are stabilized if specific external factors are included.
As observed in our study, the dynamics described by the different discrete versions is similar within an interval of confidence, but special care has to be taken for larger values of the parameters.Hence an important conclusion of our study is that a single discrete model is not enough to make assertions about a population.Thus we suggest that to check the re-ability of a continuous model is convenient to check several versions of discretizations of the model and base conclusions in the interval of confi- dence.Moreover, we may include all discretizations in the discrete map proposed in (9).

FIGURE 6 5 FIGURE 7
FIGURE 6 Bifurcation diagram for exponential model with cubic migration.