Logistic difference equation

Logistic difference equation. Contributed by Sebastian Bonhoeffer; adapted for BioSym by Stefan Schafroth. In a influential paper in 1976 the Australian theoretical ecologist Robert May showed that simple first order difference equations can have very complicated or even unpredictable dynamics. Here we explore the route into chaotic behaviour using the Logistic Difference Equation (LDE) as a model. In particular, we will address how chaotic dynamics may be characterized. 1. Introduction. The LDE can be coded in a few lines of Matlab code. This function can then be used to study the dynamic behaviour as a function of the parameter r describing the population growth rate. The dynamical behaviour can be investigated by studying time courses, phase diagrams, and bifurcation diagrams. The following concepts and methods that can be illustrated with the LDE model: Chaotic dynamics Periodic dynamics Bifurcations Definition of chaos Simulation of simple single-species difference equations Bifurcation diagrams Phase diagrams. The important message is that simple models can show complex behaviour. Hence, if we observe a complex dynamic behaviour of a natural population in the field, this does not necessarily imply that the underlying rules determining the dynamics have to be complex too. Moreover, the complex behaviour of the logistic difference equation is in strong contrast with the comparatively simple behaviour of the logistic differential equation. As a rule of thumb difference equations show more complicated behaviour than the corresponding models in continuous time. 2. Mathematical Model. 2.1 Difference equations versus differential equations. Frequently population dynamical models assume continuous time. Models in continuous time are usually appropriate for organisms that have overlapping generations. However, many biological populations are more accurately described by non-overlapping generations (such as insect populations with one generation per year or annual plants). The dynamics of these populations are often more appropriately expressed by so-called difference equations. Here, the population size of the next generation t+1 is expressed as a function of the population size in the current generation t. In general a single species difference equation can be expressed as. In analogy to the exponential equation the simplest difference equation is obtained if C ( n t ) is a constant: It is easy to see (by iteration of the above equation) that the solution is given by. where n 0 is the population size in generation t = 0. This equation is the analogue of the exponential function in discrete time. 2.2 Logistic difference equation. The logistic difference equation is given by. It can be derived as a discrete time analogy to the logistic differential equation, which is given by. where x is the population density (scaled by its carrying capacity) and r is the maximal growth rate of the population at low values of x. Another equation that is often referred to as the logistic difference equation or logistic map is given by x t + 1 = r x t ( 1 - x t ) , where 0 ≤ r ≤ 4 and 0 ≤ x ≤ 1. This equation displays analogous dynamical behaviour as Eq. 3. 2.3 The route to chaos as illustrated by the bifurcation diagram. Figure 1 shows a bifurcation diagram. This plot is obtained in the following way. For each value of r, the logistic difference equation is iterated for n steps (starting from a random initial number) to attain stable behaviour (if there is any). Then a further m iteration steps are performed and x t is plotted at the resulting time points n x t = 1 for large t (which implies that the population nt is always at its carrying capacity K). At r=2 there is the first bifurcation (period doubling) as the solution begins to oscillate between two values for this value of r. As r increases further there are further period doublings, resulting in cycles of period 4, 8, 16, and so on. Eventually, for even larger values of r the logistic differential equation shows chaotic behaviour, which means that the population behaviour cannot be predicted accurately for longer periods of time. Two trajectories starting from nearly identical values will diverge further and further away from each other. 2.4 Phase diagrams and stable cycles. Another way to explore the dynamical behaviour of the logistic difference equation are phase diagrams in which the value of x t + k is plotted against x t for a given time series and a given choice of k. Figure 2 shows the phase diagrams for several values of k. If (after a transient phase of equilibration) x t + k = x t for all values of t, then the system shows stable cycles with period k. (The bifurcation diagram 1 indicates that there are parameter values of r for which the system shows stable cycles). If for a given value of r the system converges to a stable cycle with period k, then the trajectory should converge to a point in the phase diagram plotting x t + k against x t . 3. Implementation in Matlab/Simulink. 3.1 Matlab script for simple difference equation. Download Matlab script for the simple difference equation (1). Try to modify it to implement the LDE (3). As this is the first time we deal in detail with the Matlab script language we will take a closer look at the code. Each Matlab function starts with the function declaration that defines the function name, the function parameters and the function results. Our main function has no parameters and no return values: The lines starting with a precent sign are comment lines. They describe what the function is doing and, if any, what the parameters and return values are. The next few lines represent the body of the function and describe what the function really does. The first statement calls another function called model that takes three input arguments (the numbers 1.2 , 1 and 10 ) and returns two result values: time and n . The plot command plots n versus time . The xlabel and ylabel commands define the labels of the x- and y-axis of the plot, respectively. In the code above we used a function model which must now be defined. This function implements the single species difference equation n(t+1) = R*n(t) which is also stated in the two commen lines of the function description. The function takes three input arguments: r the growth factor, n0 the initial population size and time the maximum number of time steps to be performed. The function has two return values named out_time and out_n . As always the function declaration is followed by the function body that contains the statements that describe the algorithm that is executed. We start by preparing two vectors: out_time is a vector filled with values 1,2. time and out_n is a vector the same size as out_time by initially filled with zeros. It will later hold the values of the population, n(t). We can now start the iteration by assigning n0 to an auxiliary variable n and then executing a loop that iterates time times. Within this loop, the current population size n is multiplied with the growth rate r and assigned to the variable n again. This way the variable n actually holds the new value of the population, the value which is the starting point for the next iteration. In addition, the value of n is stored in the output vector out_n at index t . 3.2 Matlab GUI program for bifurcation plots. Download these two files: These files represent a Matlab program with graphical user interface (GUI) that allows you to experiment with the bifurcation diagram (see Figure 1) of the LDE (3). You can change the r range by adjusting the sliders and/or enter the values into the edit boxes. 3.2 Simulink models. 3.2.1 Simple difference equation. Download Simulink model for the simple difference equation (1). The model is shown in Figure 4 and the result of an example run is shown in Figure 5. We introduce here a new Simulink block, the unit delay . This block outputs at time t the value that was applied at its input at time t-1. With this block we are able to implement the difference equation (1). 3.2.2 Logistic difference equation. Download Simulink model for the logistic difference equation (3). The model is shown in Figure 6 and the result of an example run is shown in Figure 7. We introduce here a couple of new Simulink blocks: Bias : This block adds a constant to the input value, in the example -1 is added. Function block : Applies a built in function on the input value. Possible functions are exp, log, cos, sin etc. Please consult the Matlab help for more details. 4. Interesting questions that you can investigate. What types of dynamical behaviour can you detect in the LDE? What defines chaotic behaviour? How can one differentiate between stochasticity and chaos? Bifurcation diagrams: Are there windows of periodic behaviour within chaotic regimes? What periodicity do they have? Advanced question: What happens if you expand the LDE to a spatial simulation (like the spatial Nicholson-Bailey model)? Is the behaviour still chaotic? 4.1 Basic exercises. Eb1 ★ : Generate bifurcation plots as in Figure 1. Amplify interesting sections of the plot by constraining the range of r and using smaller increments. Generate time series plots for interesting values of r. Eb2 ★ : Generate chaotic time series (e.g. take r = 3.5) and look at the distribution of x values attained in the series. Is it completely random (i.e. uniformly distributed)? Eb3 ★★ : Are there windows of periodic behaviour within chaotic regimes? Try to develop a method that computes the periodicity of a time series. (A periodicity of k means that x t+k = x k for all t that are sufficiently large such that the dynamic behaviour has stabilized). Are there parameter values of r for which the LDE shows dynamic behaviour with uneven periodicity? Advanced option: plot periodicity as a function of r. 4.2 Advanced/additional exercises. Ea1 ★★ : Plot the dynamic behaviour of the LDE in a phase diagram (see section 1.5 and Figure 2). When are periodic solutions stable? What are the graphical conditions for stability? (Hint: look at the slopes of the intersecting lines.) Ea2 ★★ : How can one define chaotic behaviour? What characterizes chaotic versus non-chaotic regimes? Try to develop a quantitative measure that allows you to distinguish between both regimes, and approximate the critical value of r at which chaotic behaviour appears. Ea3 ★★★ : Develop a spatial extension of the LDE (along the lines of the spatial Nicholson-Bailey model, see corresponding module). Assume that a certain fraction of individuals migrates to neighbouring fields. How does space affect chaotic dynamics? Note that this may be a work-intensive (but fun) extension to the LDE! 5. Glossary. Difference equations: Difference equations describe the dynamics of variables (i.e. the population size) in discrete time, whereas differential equations describe the dynamics in continuous time. Bifurcation diagram: a plot of the fixed points of the system as a function of a parameter. Most interesting are bifurcation points at which fixed points may be lost or born, or may lose or attain stability. Phase diagram: In the context of a single species difference equation a phase diagram displays a time series by plotting x[t+k] against x[t] for a chosen value of k. More generally a phase diagram plots a time series (trajectory) in phase space in which the dimensions are the variables of the dynamical system. Periodic behaviour: A discrete time series is periodic with period k if x[t+k] = x[t] for all t.