Constrained Population Growth

The study of state-dependent population growth rates traces to the mid-1800s. Pierre-François Verhulst famously proposed that the constraints on population growth increase as the population increases. [Verhulst-1845-NMAR] (p.8) proposes that the limited availability of good land might explain these changes. With increasing crowding, the birth rate eventually decreases, and the death rate eventually increases.

This observation readily links Verhulst’s work to modern discussions of environmental carrying capacity, which is the maximum population of a species that a particular habitat can sustain. However, Verhulst did not directly consider the concept of carrying capacity, speaking instead of a natural long-run level of the population. The present chapter adopts a different terminology, the sustainable population, which later sections connect to the mathematical concept of a steady state.

Logistic-Growth Map

The corresponding evolution function is a discrete-time Verhulst-Pearl growth function. For any given \(\bar{g}\), this is more commonly called a logistic-growth function or logistic-growth map. [1] Consider the following binary real function. For any given \(\bar{g}\), it transforms a real number (the input) into a real number (the output). Mathematically, this function certainly can be defined on all of the real numbers. However, the conceptual model makes it possible to be more specific: the concept of relative population only makes sense of nonnegative real numbers.

Model Parameters under Logistic Growth

Recall that in simulation modeling, a model parameter is a constant in the model implementation that can affect the simulation outcomes. In the PopulationGrowth01 project of the Exponential Growth chapter, the population projections depend on two parameters: the constant population growth rate, and the initial population. These are the two model parameters in the project.

In the LogisticGrowth01 project of the present chapter, the initial population is again a model parameter. The growth rate is no longer constant, but a comparable role is now played by the maximum growth rate (\(\bar{g}\)) parameter. In addition, the logistic-growth model adds third parameter, reflecting a new consideration. Population growth now takes place in an environment characterized by a maximum sustainable population (\(P_{ss}\)). This planetary carrying capacity is a new model parameter.

Population Projections

Recall from the Exponential Growth chapter that following the setup phase, a typical discrete-time simulation includes an iteration phase. The iteration phase repeatedly executes the simulation schedule. At the core of this schedule is the model update, which implements the evolution rule for the simulation. (This is the computational implementation of the evolution rule derived from the conceptual model.) Each iteration of the model schedule updates the model state.

In the current project, the model state is simply the relative population. The model update uses the nextVP function to update the relative population (code:relpop). Nevertheless, the ultimate goal is generation of a sequence of population projections. This constitutes the next component of the LogisticGrowth01 project: producing the population projections. As with the exponential projections, the logistic population projections will naturally depend on the initial population. However, as shown below, population growth now follows a logistic trajectory. As a result, that dependence diminishes over time.

Approach this population simulation by implementing the following trajectory function sketch. Here, the specification () -> Sequence indicates that the function has no input parameters but returns a sequence. This is the trajectory of population projections. Note that in this sketch, the function is not pure: this sequence it produces depends on its environment. Specifically, it depends on the values of the model parameters. These parameters occur in the proposed trajectory function as free variables; their values come from outside the function definition.

function signature:

trajectory: () -> Sequence

summary:

• Do any needed setup of the simulation model.

• Compute the initial relative population (relpop=P0/Pss).

• Repeatedly (\(50\) times) use nextVP to update the relative population, making use of the parameterized maximum growth rate.

• Transform each relative population into an absolute population (relpop * Pss).

• Return a sequence containing the entire trajectory of population projections.

Nonlocal Variables

On the one hand, relying on pure functions is often a good practice: they are easy to understand and have completely predictable behavior. On the other hand, it is desirable for function definitions to avoid hard coding values that we may need to change. The present chapter adopts one possible resolution of this tension: function definitions sometimes may usefully refer to nonlocal variables. (Later chapters explore other approaches.)

A function definition depends on a nonlocal variable when there is an outside influence on the values returned by the function. The value of a nonlocal variable exists outside of the function definition. This implies that the function behavior depends on a potentially changeable part of the environment in which it is executed.

A pure function completely characterizes the relationship between its inputs and its outputs, with no outside influences. In contrast, when a function depends on a nonlocal variable, the relationship between a function’s explicit inputs and its output is unknown until the value of the nonlocal variable is specified. Pure functions are easier to completely understand and therefore are safer to use. Nevertheless, at times impure functions prove convenient, especially early in the development of a program. The rest of this section illustrates one possible use of nonlocal variables in the conduct of simulation experiments.

Working with the Trajectory

At this point, we have developed all the components needed to run the simulation and collect data from it. The setup phase of the simulation parameterizes the computational model and performs any other needed initialization. It may additionally prepare to collect simulation data. The iteration phase generates a population trajectory. It may additionally collect simulation data. It is natural to implement an activity that handles both phases.

Population-Projection Plots

Once we can access the population projections, it becomes natural to plot them. Plots often make it easy to see patterns in data that might be challenging to detect in raw numbers. This constitutes the next activity. The following plotTrajectory activity summary leaves out many details, since the core of this activity is simply to produce a plot.

activity:

plotTrajectory

summary:

Produce a run-sequence plot of the population trajectory produced by your LogisticGrowth01 project.

Here are some possible ways to execute the plotTrajectory activity. One may use spreadsheet software to import the exported data and then plot it. One may use other plotting software to plot this data. Or, one may use builtin plotting facilities of the programming language or toolkit used for the simulation. This last approach can use the simulated sequence directly, rather than importing the data.

Execute the plotTrajectory activity for the LogisticGrowth01 project. Compare the new population plot to the plot from the PopulationGrowth01 project.

Figure f:CompareExpLog presents a line chart based on the LogisticGrowth01 project along with the one from the Exponential Growth chapter. The difference between the projection based on exponential growth and the projection based on logistic growth is immediately evident. In the logistic-growth model, the population growth rate slows over time, while in exponential-growth model it does not. Furthermore, over time the population in the logistic-growth model appears to be constrained by the fixed sustainable population.

Role of Maximum Growth Rate

The previous section developed a logistic-growth simulation model and produced a set of population projections. These projections eventually diverge substantially from the exponential growth projections of the Exponential Growth chapter. The 50-year forecast from the logistic-growth population model is much lower than the exponential-growth forecast. The new forecast appears more realistic. For example, it better matches the median-variant projection of the United Nations.

The present section considers a natural question about the logistic-growth model: how sensitive are these population projections the model parameterization? This constitutes a research question that may be addressed by means of a simulation experiment. In particular, the maximum population growth rate (\(\bar{g}\)) appears to be a key model parameter for these projections. This section provides a first approach to the question of how changes in the maximum-growth-rate parameter alter the population projections. The approach is experimental: we simulate the population trajectories for various maximum growth rates and then compare them.

Primitive Experimentation

A model parameterization, or scenario, gives a concrete value to each model parameter. During a single simulation run, all of the model parameters retain the same fixed values. In a deterministic system-dynamics model without exogenous inputs, a specifying a model parameterization fully determines the outputs of the simulation.

A simulation experiment compares multiple scenarios. This means that some parameters vary across simulation runs. Simulation experiments allow exploration of a simulation model by exposing how parameter variations affect the system behavior. The parameters selected for variation depend on the particular questions that motivate the experiments. These are the focal parameters of the experiment.

The maximum growth rate is the focal parameter for the present experiment: we wish to discover how that value of \(\bar{g}\) affects the simulation outcomes. Unfortunately, in our LogisticGrowth01 project, this value is hardcoded in the source code. Nevertheless, as a first approach to experimentation in this model, briefly try out the following error prone and tedious strategy. Repeatedly edit the project code to change the value of \(\bar{g}\), rerunning the simulation each time. After each run, recreate the associated run-sequence plot.

The point of this exercise is to gain insight into the limits of this kind of manual approach to simulation experiments. So as you explore this experimentation strategy, you may want to write down some notes about any advantages or disadvantages that you discover.

Hard Coded vs Named Parameters

In the Exponential Growth chapter, the definition of the exponential growth function included a hard coded 1% population growth rate. This ensured that nextExponential was a pure function; the function definition provides all the information needed to calculate the function’s output for any input. Functional purity is a very handy: for any suitable input, there is never any confusion about what this function should produce as output (subject to the limits of the computer). The binary nextVP function of the present chapter is also pure: given two real inputs, it deterministically produces a real output. It is closed to all outside influence, and it has no side effects.

However, the LogisticGrowth01 project hard codes the model parameters. As illustrated by the previous exercise, for the current simulation experiment, this has drawbacks. Although the maximum growth rate must remain constant during any simulation run, the experiment requires it to change across simulation runs. A hardcoded model parameter becomes an annoyance whenever a simulation experiment requires repeatedly resetting it to varied values. It is wise to assume that any attempt to coordinate such code changes eventually will fail. To address this, let us refactor the trajectory producing function, making it less dependent on its execution environment.

Refactoring the Trajectory Computation

Specifically, use the following function sketch to implement a trajLG02 function, which handles the trajectory computation. This is a ternary function: it accepts the three model parameters as inputs, and produces a simulated logistic-growth population trajectory. (The function summary provides descriptive guidance, but produce the population trajectory however you wish.)

function signature:

trajLG02: (Real, Real, Real) -> Sequence

parameters:

• \(\bar{g}\), the VP maximum growth rate parameter

• \(p_\text{ss}\), the sustainable population

• \(p_0\), the initial population

summary:

• Compute the initial relative population, \(r_0 = p_0/p_\text{ss}\).

• Partially apply \(\mathrm{nextVP}\) to produce \(f = r \mapsto r + \bar{g} \, r \, (1 - r)\).

• Iteratively (\(n\) times) apply \(f\) to produce the relative-population trajectory \(\langle r_0, f[r_0], \dots, f^{\circ 50}[r_0] \rangle\).

• Return the corresponding population trajectory. (Each population scales a relative population by the sustainable population: \(r * p_\text{ss}\).)

When computing a population trajectory, the trajLG02 function depends on the value assigned to \(\bar{g}\). Applying this function to a new value of \(\bar{g}\) produces a different population trajectory. Each change in \(\bar{g}\) produces a new trajectory. The trajLG02 function provides an abstract representation of all of these trajectories. This allows us to easily experiment with changes in the maximum growth rate. Take advantage of this by executing the following collectTrajectories activity.

activity:

collectTrajectories

summary:

Produce a useful collection of three logistic-growth population trajectories, one each for the following values of \(\bar{g}\): 1%, 3%, and 6%. Use trajLG02 to produce each trajectory. Keep the sustainable population and initial population at the baseline values, given in Table t:LogisticParameters.

Abstractly, producing the collection of trajectories is effectively a mapping operation, using a partially applied trajLG02 function over a sequence containing the maximum growth rates. The result of such a mapping operation is a sequence of trajectories. However, you may store the three trajectories in any way you find most helpful for the following exportTrajectories activity. which is to export the trajectories to a CSV file.

activity:

exportTrajectories

summary:

After executing the collectTrajectories activity, write the three population trajectories to a columnar CSV file named LogisticGrowth02.csv. That is, put one trajectory in each of its three columns.

Finally, consider the following plotTrajectories activity. This activity can also benefit from the collectTrajectories activity. After producing the collection, use the lessons learned in the LogisticGrowth01 project to plot each trajectory.

activity:

plotTrajectories

summary:

After executing the collectTrajectories activity, produce a run-sequence plot for each trajectory. If possible, place the plots in a single chart, along with a helpful legend.

Experimental Results

Figure f:varyMaxPopGrowth illustrates typical results from completing this exercise. Raising the maximum growth rate produces initially more rapid growth, so the corresponding population projections are higher. Fifty-year population projections are clearly sensitive to the maximum growth rate in the logistic model of population growth.

Steady State

The LogisticGrowth02 project demonstrated that, over the chosen forecast horizon of 50 years, fairly small changes in the growth rate parameter have fairly large effects on the population projections. In this sense sense, population projections are sensitive to this parameter. However, longer term projections eventually converge to the sustainable population. In this sense, the population projections are eventually similar, despite these changes in the parameterization. Although changing the maximum growth rate changes the speed of convergence, it does not change the sustainable population. This implies that in the LogisticGrowth02 project, a long-run population projection does not depend on the maximum growth rate.

This observation leads us toward the concept of a steady state of a dynamical system. A steady state is a state that, in accordance with the evolution rule of the system, does not change over time. If a deterministic dynamical system is in a steady state, it stays there.

In LogisticGrowth02 project, the sustainable population (\(P_{ss}\)) determined the long-run outcomes of the logistic-growth simulations. Furthermore, \(P_{ss}\) is a steady-state population. That is, if the population equals \(P_{ss}\), it will not change. Restated in terms of populations relative to the sustainable level, if the relative population equals \(1.0\), it will not change.

Steady States as Fixed-Points

From a system-dynamics point of view, it is natural to speak of a steady state of a system. The related concept from mathematics is a fixed point of a function. When applying a function to some value reproduces that value, the value is a fixed point of the function. So a steady state is a fixed point of an evolution function.

Recall our logistic-growth function for the relative population: \(p \mapsto p + g \, p \, (1 - p)\). Iterating this function produces a trajectory for the relative population. It is easy to show that this system has a steady state, where the output value of this function is the same as the input value. The steady-state values solve the following equation in \(p\).

This equation expresses the requirement that the function produce a value equal to its input. Given \(g \neq 0\), we can see by substitution that both \(0.0\) and \(1.0\) are steady-state values of this dynamical system. That is, they are fixed points of the logistic-growth function.

Recall that a relative population of \(1\) corresponds to a population that equals the sustainable population. A system that attains the sustainable population stays there. As we might guess intuitively, a population of \(0\) is also a steady state: an extinct species stays extinct.

The LogisticGrowth02 project produces \(50\) years of population projections for various maximum growth rates. In each case, longer projections eventually converge to a relative population of \(1\). None of these trajectories converge to 0. This convergence behavior indicates that these two fixed points have very different dynamical properties.

In the LogisticGrowth02 project, the value \(1\) is a stable fixed point of the logistic functions: when the population is near the sustainable population, it approaches it. Equivalently, the value 1 is an attracting fixed point, or attractor, in these dynamical systems. In contrast, the value 0 is an unstable fixed point: when the population is a bit above \(0\), it continues to grow. Equivalently, we say that the value 0 is a repelling fixed point, or an repellor, of these dynamical systems.

A New Experiment

Despite its simplicity, the Verhulst-Pearl logistic growth function can produce more complicated behavior than monotone convergence toward the steady state. As a first observation to support this claim, consider the possibility of overshooting the steady state when \(g > 1\). For example, confirm that \(\mathrm{nextVP}[1.5,0.8]=1.04\), exceeding the steady-state value of \(1.0\). A subsequent iteration then produces a value below the steady state. It seems that at large values of \(g\), a logistic-growth function can produce a trajectory that fluctuates around a steady state. This is the case, as demonstrated by the following exercise.

To support this exercise, and more generally to allow varying the number of iterations, slightly refactor the trajectory function to take the following three arguments: the maximum growth rate, an initial value, and a number of iterations. Call this new function iterVP.

function signature:

iterVP: (Real, Real, Integer) -> Sequence

parameters:

• \(g\), the VP growth parameter

• \(p_0\), the initial population

• \(n\), the number of iterations

summary:

• Let \(f = p \mapsto p + g \, p \, (1 - p)\).

• Iteratively (\(n\) times) apply \(f\) to produce the trajectory \(\langle p_0, f[p_0], \dots, f^{\circ n}[p_0] \rangle\).

• Return the entire trajectory.

Experimental Results

Figure f:varyMaxPopGrowth2 illustrates typical results from completing this exercise. The now larger values of the maximum growth rate produce quite different behavior than before. Letting \(g = 1.9\) still produces a convergent trajectory, but it fluctuates around the steady state. Raising \(g\) to \(2.0\) still produces a convergent trajectory, and it still fluctuates around the steady state, but it converges more slowly. (Much more slowly, as it turns out.) Raising \(g\) still further to \(2.1\) produces a new behavior. The trajectory still fluctuates around the steady state, but it does not appear to converge.

Picturing Fixed Points

Why does that last case (\(g = 2.1\)) appear not to converge? Let us explore this with a common visualization of the fixed points of this logistic function. Create an ordinary function plot in the Cartesian plane, along with a \(45^\circ\) line. At its fixed points, a function plot crosses the \(45^\circ\) line. In this case, the function is the evolution function of a dynamic system, so these points of intersection identify the possible steady states. Figure f:logistic-vp210 provides an illustration.

As expected from the algebraic discussion of the fixed points of the nextVP function, we always see the two lines intersect at \((0,0)\) and \((1,1)\). This remains true even at much higher values of gmax; the locations of the fixed points do not depend on gmax. However, this does not mean that the behavior of the model is unchanged. Now each fixed point is unstable, because the function slope is too steep.

It is natural to be a bit puzzled that the steady state is no longer an attractor. This results from the steeper slope of the logistic-growth function. For a fixed point to be an attractor, the slope of the function at that point must be less than one in absolute value. Increases in \(gmax\) increase the slope of the logistic-growth function at its fixed points. The sustainable population is still a steady state for the system, but it is no longer a stable steady state. To be a stable fixed point, the slope must lie in \((-1\,..\,1)\), but now the slope is \(3.25\) at \(p=0.0\) and \(-1.25\) at \(p=1.0\). Raising the maximum growth rate makes the function steeper and steeper at the two steady states, until eventually the steady state at \(1.0\) is no longer stable.

Periodic Attractors as Fixed Points

So the steeper slope explains the new lack of convergence, but what explains the apparently periodic behavior? To approach this question, suppose this function actually does have periodic values, recurring every two periods. Begin a function iteration with one of these points, apply the logistic-growth function to it, and then apply the logistic-growth function to the result. This reproduces the original point. That is, such a periodic point is a fixed point of the twice iterated function.

As a convenient notation, for any unary function \(f\), represent the twice iterated function as \(f^{\circ 2}\). Now let \(f = p \mapsto \mathrm{nextVP}[2.25, p]\), so that \(f^{\circ 2}\) is a twice iterated logistic function. We may then use the same graphical approach as before to reveal the fixed points of this iterated function. Simply plot \(f^{\circ 2}\), and search for its intersections with the \(45^\circ\) line.

Of course any fixed point of \(f\) must also be a fixed point of \(f^{\circ 2}\). So the new plot must rediscover the original steady-state points. And as shown in Figure :name: f:logistic-vp2x225, the iterated function \(f^{\circ 2}\) indeed has fixed points at \(0\) and \(1\). But in addition, there are now \(2\) more intersections: the two periodic points (near \(0.71\) and \(1.17\)).

This is an example of a periodic attractor, and the two additional fixed points are called periodic points. More specifically, the two points can alternate in a period-2 cycle. What is more, this two cycle is attractive. We see that visually because the slope of \(f^{\circ 2}\) at the periodic points is clearly less than unity. So if we start at some other value of \(p\), the trajectory will approach the period-2 cycle over time. This is an example of a limit cycle: the cycle is the limiting behavior of the trajectory.

Longer Cycles

Our discovery of period-2 cycles shows that our iteration of our simple logistic-growth function can generate surprisingly complicated dynamics. It is natural to wonder if longer cycles are possible. For example, does \(f^{\circ 4}\) have any additional fixed points? Letting \(g = 2.25\), a quick plot reveals no additional fixed point. Nevertheless, there is more to come.

Another Experimental Result

Figure f:logistic-gmax250 illustrates typical results from completing this exercise. Letting \(g = 2.5\) produces a fluctuating trajectory around the nonzero fixed point, but this time there appears to be a 4-cycle.

Letting \(f=p\mapsto nextVP[2.5,p]\), one may try to shed light on this with a function plot. The plots of \(f\) and \(f^{\circ 2}\) are easy to interpret: the 1-cycles and 2-cycles are clearly unstable (due to the steep slopes there). Making a useful plot of \(f^{\circ 4}\) is challenging: While it is possible to make out the four points of the 4-cycle, this is not easy unless the plot is very big (or zooms in). It is time for a better approach.

Period Doubling

We found that when the VP growth parameter is \(2.10\), the trajectory is attracted to a period-2 limit cycle. Two points have long-run relevance for this system, which alternates between them. When we increase the VP growth parameter to \(2.5\), the trajectory is attracted to a period-4 limit cycle. Four points have long-run relevance for this system. In the face of this discovery, it is natural to ask whether such period-doubling continues to result from increases in the VP growth parameter. The answer is yes—and with faster and faster period doubling. (This is called a period-doubling cascade.) When we raise gmax to \(2.55\), the trajectory is attracted to a period-8 limit cycle, but unfortunately this is almost impossible to discern with such plots. It becomes harder and harder to detect these oscillations visually.

Brute Force Computation of Limit-Cycle Points

In addition to such basic problem of visual discrimination, regularities in the limiting behavior of a function iteration may be concealed by the early behavior. To deal with this problem, simulation analysis often incorporates a burn-in period. This involves discarding initial iterations of a dynamic system in order to better focus on the behavior it ultimately settles into. The burn-in period is the number of discarded iterations.

There is another problem: computers generally work with approximate numbers. If computers were capable of infinitely precise computations, the actual trajectory would approach its limit cycle over time. Unless the initial point is a periodic point, however, the process of convergence would continue forever. However, given the limitations of digital computers, an actual cycle between approximate numbers emerges after a finite number of iterations. This suggests a strategy for numerically finding periodic points: run the simulation for a while, and then search for repeating values towards the end of the trajectory. That is, discard the first part (“prefix”) of a long trajectory, and then examine the final segment for periodic behavior. One approach to producing such a post-burn-in sequence segment is illustrated by the following segmentVP function sketch.

function signature:

segmentVP: (Real, Real, Integer, Integer) -> Sequence

parameters:

• \(g\), the VP growth parameter

• \(p_0\), initial value for trajectory

• \(n_1\), the length of the trajectory prefix to discard

• \(n_2\), the length of the trajectory segment to keep

summary:

• Let \(f = p \mapsto p + g \, p \, (1 - p)\).

• Starting with \(p_0\), iteratively (\(n_1\) times) apply \(f\) to produce the value \(p_{n_1} = f^{\circ n_1}[p_0]\).

• Return iterVP[g,p_{n_1},n_2 - 1].

This is a good start, but if this detects a limit cycle, much of the sequence segment is likely to be redundant. The limit-cycle points are the set (i.e., no repetitions) of the elements in the sequence segment. In principle, discarding the duplicates is fairly trivial, since most language include some facility for easily doing this. In practice, however, we again run into the limitations of floating-point computations: multiple nearly equal approximate numbers may represent a single limit-cycle point. For this reason, let us treat numbers that are very similar as actually being the same. One way to do this is to discard the final decimal places of each computer number prior to comparing them. The following sketch of a nudupes function captures this idea.

function signature:

nodupes: RealSequence -> RealSequence

parameters:

• \(\mathrm{xs}\), the input sequence of real numbers

summary:

• Let \(\mathrm{ys}\) be the values of the \(\mathrm{xs}\), but rounded to 12 decimal places.

• Return the unique elements of the \(\mathrm{ys}\).

A bifurcation is the doubling of the period of the limit cycle. A bifurcation plot is popular way to illustrate the effects of changing the parameter gmax. For each value of g, it plots the limit-cycle points determined as above (but often with longer trajectories). It displays the points that are occur in the limit cycle. This will be only 1 point for \(g < 2\), but will then move into a period-doubling cascade, including more and more points until a visual display of all the points becomes impossible.

Once we choose the values of the VP growth parameter to consider, we are ready to create a bifurcation plot. This is a scatter plot showing the periodic (limit-cycle) points for each value of the VP growth parameter. This chapter describe a brute-force approach to producing this plot, while acknowledging that this has a number of problems. For example, how do we know if we have chosen an adequate burn-in period? And, how do we know if we have chosen an adequate segment length? Nevertheless, brute-force approaches to the creation of a bifurcation diagram remain far the most common. Experimentation indicates that the following brute-force approach should quickly produce a revealing plot for this logistic-growth model. Figure f:periodicPointsLogistic illustrates a typical result of this brute-force approach.

On to Chaos

Examining the bifurcation diagram, we can easily detect the first three period doublings discussed above. But then the situation quickly becomes much messier. In fact, above about \(2.57\), many situations arise where there is no periodicity at all. This is the realm of chaotic dynamics. The behavior of this simple nonlinear dynamical system is no longer just complicated, it is complex.

This complex behavior has spawned a very large literature. Of particular interest is that, when dynamics are chaotic, this simple, deterministic, nonlinear system becomes effectively unpredictable. To pin this down a bit, it means that even miniscule changes in the initial conditions can lead to utterly divergent long-run behavior. To put it more provocatively, long-run prediction effectively becomes infeasible.

Conclusion

This book repeatedly considers four steps: developing a conceptual model, implementing it as a computational model, using the computational model to run a simulation, and examining the data produced by the simulation. This chapter illustrated these components of simulation modeling by developing a simple model of logistic population growth.

Supplementary Exploration

1. Use a spreadsheet to produce a 50-year population projection, under the same population assumptions as above. In the Spreadsheet Introduction supplement, find details for this spreadsheet exercise in the discussion of step-by-step population projection. Compare your spreadsheet results to the simulation results above.

2. Change the trajectory function in the LogisticGrowth01 project so that is it a pure function, which consumes three arguments. (You may keep the hard coded trajectory length.) Make the implied changes elsewhere in the project code. (Each use of trajectory now needs three arguments.)

3. Each step, in addition to writing the current population to LogisticGrowth01.csv, write the current relative population. Write the two numbers on a single row, separated by a comma. (This data is in CSV format.)

4. Reproduce Figure f:compareExpLog, which compares the outcomes under exponential growth and logistic growth.

5. Reproduce Figure f:varyMaxPopGrowth, which compares the logistic-model outcomes under three different values of gmax. Predict what will happen if you extend the simulation another \(50\) years: will the final predictions be even further apart, or will they be closer together? Explain the reasoning behind your hypothesis. Run simulations to test your hypothesis.

6. Impose the constraint that the input to the nextVP function is nonnegative. (The best approach to this depends on the programming language.)

7. Create a slider for the initial population (\(p_0\)). Use this slider for simple model exploration. Does the initial population matter for the long-run behavior of the system?

8. Instead of hard coding the output file name in the exporting activity, introduce it as a global constant (say, outputFile).

9. (Advanced) In the LogisticGrowth02 project, instead of hard coding the number of iterations in the trajLG02 function, make it another function parameter (say, nIter). How does this affect how you produce your plots?

10. (Advanced) In the Spreadsheet Introduction supplement, do the regression exercise on estimating the sustainable population for the U.S.

Additional Resources

Global Variables

Recall that simulation models in this book typically include a setup phase that sets the model parameters and initializes any global variables. The initial values of the global variables often depend on the model parameters.

At this point, there are two global variables in the computational model. The first is relpop, which holds the current population relative to the sustainable population. (So relpop is the source-code variable that corresponds to the \(p\) in the mathematical discussion above.) Given the values of model parameters, initialize relpop to P0 / Pss. Second, as a convenience, the model includes a global variable named ticks. This variable should track the number of iterations of the model step, so initialize ticks to \(0\) during the setup phase.

In addition to setting the values of the model parameters, the setup phase includes the initialization of any global variables. For example, the exponential growth model of the Exponential Growth chapter introduced the current population as a global variable of the model. This offers a simple way to track the state of the population during a simulation. In the present chapter, a similar role is played by the relative population (relpop). This is the ratio of the population to the sustainable population. Whenever it might be needed, the actual population may be recovered as Pss * relpop.

Incremental Data Export

An alternative to end-of-simulation bulk export is incremental export of the data. This is particularly advantageous when a simulation generates very large amounts of data. For incremental export, the setup phase includes preparing an output file and writing the initial data. The iteration phase then appends data to this output file after each simulation step.

For incremental export, the setup phase typically includes a little file preparation. For example, the setup phase typically discards any existing data from prior runs of the simulation. In addition, for incremental data export, the setup phase typically writes the initial model state to the output file. The following summary of a setupFile activity captures these considerations.

activity:

setupFile

summary:

• Set up an empty output file.

• Write the initial population as plain text to the first line of the empty output file.

For incremental export, the end of the iteration phase should write new data to the output file. Even though the data-export step is currently very simple, this section introduces a recordStep activity to handle it. For the current model, execution of the recordStep activity appends a number to the output file by writing it on its own line. That number is the current population.

activity:

recordStep

summary:

• Open the output file for appending.

• Append the current population to the file, on its own line.

• Close the output file.

Naming the Output File

The earlier discussions of file handling assumed that the name of the output file is hardcoded. This is fine when a deterministic simulation runs under a single parameterization. Multiple simulation runs create some new considerations when it comes to data export, since each simulation run will overwrite the last output.

For example, suppose each simulation writes a population trajectory to file. It would be possible after each simulation run to copy the output file to a new location, but this is again tedious and error prone. Two possible responses are either to alter the code to run the entire experiment and write all the resulting projections to a single file or to ensure that each parameterization produces a unique output file. For the moment, consider the latter.

When a parameter changes across simulation runs, it should be possible to somehow link the parameterization to the output file. This book explores a number of approaches to this problem, but here is the simplest: specify the parameterization in the name of the output file. This eliminates one source of potential confusion in the analysis of the results of the initial simulation experiment. Any change in the value of the gmax parameter should automatically change the name of the output file.

Revised Plot Code

Revised Slider

A goal of this section is to consider a more extensive set of values for gmax.

Revised Buttons

The runsimG activity permits very simple command-line management of the simulation, and calling it with the Run Sim button makes this even simpler for interactive exploration. Conventionally, however, two separate buttons control the setup phase and the iteration phase. So add a Set Up and Go buttons.

The Set Up button conventionally initializes the simulation model and typically sets up the GUI. Similarly, in addition to running the model’s core simulation schedule, the Go button typically ensures updates to the GUI. These details are hidden from new users of the model, who can rely on the standard GUI interface.

Visualizing Logistic Population Growth

GUI Experimentation

This section emphasizes interactive visual exploration of the logistic-growth model. It therefore focuses on the development of a useful graphical user interface (GUI). After completing this section, you will be able to create GUI buttons for simulation control. You will also be able to use sliders for interactive experimentation. The focal output of this section is a time-series plot, which provides a helpful visualization of the simulation data. A core goal is increased familiarity with the basic GUI facilities of your chosen language or simulation toolkit.

Invariance of the Fixed Points to the Growth Rate

It may be a difficult to distinguish the \(45^\circ\) line from the function plot when gmax is \(0.03\). As you increase the value of gmax, the gap between the two plots becomes more visible. Now, consider a substantially higher value of gmax, which makes it much easier to detect the two fixed points. Use the slider to set gmax to 1.6. Figure f:vpfp160 provides an illustration of the two fixed points of the logistic-growth function when gmax is \(1.6\).

New Short-Run Behavior at Higher Growth Rates

Figure f:vpts160 shows that the time series produced at such a high value of gmax: is qualitatively different: the simulation no longer produces a monotone sequence. However, increases in the number of iterations show that the population still converges to \(Pss\). The sustainable population is still a stable steady state (attractive fixed point).

References

LogisticGrowth01 (NetLogo)

Based on your reading in the Introduction to NetLogo supplement, create a new model file named LogisticGrowth01.nlogo. The Code tab of this file should hold all the code for the LogisticGrowth01 project.

to test-nextVP
  type "Begin test of `nextVP` ... "
  if (0 != nextVP 0.03 0) [error "bad next relpop fr 0"]
  if (0 != nextVP 0.01 0) [error "bad next relpop fr 0"]
  if (1 != nextVP 0.03 1) [error "bad next relpop fr 1"]
  if (1 != nextVP 0.00 1) [error "bad next relpop fr 1"]
  print "`nextVP` passed."
end

LogisticGrowth02 (NetLogo)

Begin the LogisticGrowth02 project by copying the LogisticGrowth01 project. Refactor the trajectory computation, as described in The Logistic Growth chapter. Implement the new computation as a trajLG02 function. Remember to test your work. A properly implemented trajLG02 function, which takes three arguments, must pass the following test.

to test-trajLG02
  type "Begin test of `trajLG02` ... "
  let _traj (trajLG02 gmax pss pop0)
  if (pop0 != first _traj) [error "bad trajectory"]
  if (51 != length _traj) [error "bad trajectory"]
  set _traj (trajLG02 0.0 pss pop0)
  if (pop0 != last _traj) [error "bad trajectory"]
  print "`trajLG02` passed."
end

Once your trajLG02 passes this test and you are satisfied with your implementation, compare it to the following NetLogo code. This is one possible NetLogo implementation of trajLG02.

With trajLG02 in hand, it is possible to approach the collectTrajectories activity. There are many possible approaches to this activity. In every case, we must produce the three needed trajectories, one for each value of \(\bar{g}\). Abstractly, this effectively maps a partially applied trajLG02 function over a sequence containing the maximum growth rates. The result of this operation is a sequence of trajectories. Concretely, even though map is a builtin function, a new NetLogo programmer may be more likely to build up a list of trajectories inside a for statement. This is a fine approach, but let us use map. This is usually the most idiomatic way to map a function across a list of values. The following gmaxExperiment function takes this approach, but it then goes further. Instead of simply returning the list of trajectories, it uses them to produce a table mapping experiment names to trajectories.

The gmaxExperiment function takes one argument, a sequence of maximum growth rates. It produces a trajectory for each one of these, using NetLogo’s builtin map function and this project’s trajLG02 function. Note how NetLogo’s arrow notation makes it very easy to create the function needed by map. The resulting unary anonymous function is a partially applied version of trajLG02, where the maximum growth rate is the only needed input.

Using table:from-list, it is possible to produce a table directly from a list of key-value pairs. The trajectories will be the values. In order to have useful keys, we produce a list of descriptive names, one for each experiment. Another use of map then allows us to pair each name with its trajectory, giving us the needed list of key-value pairs. We report the table produced from this list of key-value pairs.

To complete the exportTrajectories activity, we turn to the to-file command in NetLogo’s csv extension. He we support the activity with an exportTrajectories procedure, which accepts the table produced by gmaxExperiment as input. For now, to keep things simple, the output file is hard coded. (How might you produce a more informative filename?) To get the list of trajectory names, apply the table:keys function. To get the list of trajectories, apply the table:values function. To use csv:to-file to produce a columnar CSV file, we first need to transpose the list of trajectories. Here we use the unzip function from the fp extension, as discussed in the NetLogo Programming supplement. Prepend the list of names to the new list of rows, so that there will be column headers in the CSV file. Then write the result to your file.

Finally, proceed to the plotTrajectories activity. As usual, add a plot widget to the Interface tab. However, this time in the widget dialog, delete the default plot pen and check the legend checkbox. Support this activity with the following plotLG02 procedure, which consumes the table of input trajectories produced by gmaxExperiment. (Plot setup remains unchanged.) Since there are multiple trajectories, we will use multiple plot pens—one for each trajectory. In is convenient to use the trajectory name for the pen name. Use NetLogo’s set-plot-pen-color command to set each pen to a random color, picked as one-of basecolors. Finally, plot each trajectory just as in the LogisticGrowth01 project.

It is now possible to execute the plotTrajectories activity by entering the following code at the NetLogo command line. As a refinement, add a button to the Interface tab that runs this code.

plotLG02 gmaxExperiment [0.01 0.03 0.06]

LogisticGrowth03 (NetLogo)

The LogisticGrowth03 project needs access to the trajectories produced by iteration of of the Verhulst-Pearl logistic-growth function. Provide these by implementing the iterVP function described in the Logistic Growth chapter. (This task reuses skills developed in implementing trajLG02.) Here is one approach. After partially applying nextVP to the VP growth parameter, the following implementation uses fp:iterate to do the function iteration.

Exporting the data and using other software to plot it is always an option. However, the following discussion focuses on plotting in the Interface tab. In the same way as in previous projects, add a plot widget to this project’s Interface tab. It may be named Trajectory. (Do not delete the default plot pen.) Since the new trajectories are more volatile than in the LogisticGrowth02 project, plotting multiple trajectories together in a single plot frame becomes too messy. Therefore, plot one trajectory at a time.

In order to support the plotting activity, write a plotTrajLG03 procedure that consumes the VP growth parameter (g) and produces a run sequence plot in the Trajectory plot widget. For now, plot the trajectory produced by iterVP[g, 0.75, 25]. Include whatever setup you find useful. Here is one approach.

For ease of experimentation, as described in the Introduction to NetLogo supplement, add a slider for the VP growth parameter. Name it gmax, and allow \(300\) values from \(0.01\) to 3.00. Since the gmax parameter is now an interface global variable, ensure that it is not additionally declared in globals. (Recall the discussion of interface globals in the Introduction to NetLogo supplement.) Now you should be able to use the slider to set gmax and then run plotTrajLG03 gmax from the command line to produce the corresponding run-sequence plot.

For even greater ease of experimentation, add a Run Simulation button to the Interface tab. Let this button execute plotTrajLG03 gmax. Using the new slider and button, conduct the three-trajectory experiment of the LogisticGrowth03 project. That is, one at a time, set the gmax slider to \(1.9\), \(2.0\), and \(2.1\) and the display the resulting trajectory. After noting the differences between these plots, try one more: set the slider to \(2.5\), and once again plot the resulting trajectory.

The next goal is to create a bifurcation plot. To begin this project, first determine the possible limit-cycle points for any given value of the VP growth parameter. As suggested in the Logistic Growth chapter, approach this by creating a segmentVP function and a nodupes function. Follow the chapter’s function sketches. Here is one approach.

In this example code, segmentVP relies on the fp extension, which as always must be declared at the top of the Code tab (using the extensions keyword). The first step is to partially apply the nextVP function to the VP growth parameter in order to produce a unary real function. To cheaply drop the trajectory prefix, this example uses fp:iterate-last. (This produces only that final value of the iteration.) To produce the needed segment, however, we then need to collect subsequent terms. The fp:iterate function accomplishes this exactly as in the previous iterVP implementation.

NetLogo’s remove-duplicates reporter makes it simple to produce a list of unique values. However, we wish to discard misleading noisy differences between numbers that result only from the finite precision of floating-point computations. Given a list of floating-point numbers, this nodupes function uses NetLogo’s precision function to round them to the \(12\)-th decimal place. Two values that are equal at this level of precision are treated as actually equal, so nodupes retains only one value for the two of them.

As described in the Logistic Growth chapter, the segmentVP and nodupes functions nicely support the creation of a bifurcation plot. For the production of this plot, add a second plot widget to the Interface tab of your LogisticGrowth03 project. Set the name to Bifurcation Plot, and edit the mode of the default pen to be Point. (This is the pen mode needed for a scatter plot. You may alternatively set this during any plot setup with the command set-plot-pen-mode 2.)

You can put the actual plot setup code and plotting code either directly in the plot widget or in a new plotBifurcations procedure. The discussion here focuses on the content of such a procedure. Since there are now two plot widgets, setup should use NetLogo’s set-current-plot command to activate the desired widget. It is also a good idea to execute the clear-plot command, so as to start from scratch.

Lastly, this plotting procedure should produce and plot all the desired points. Since we have no other use for these points, there is no need to collect them all (e.g., in a list). Instead, it simply plot each point as it is produced. For each considered value of the VP growth parameter (g), produce the periodic points as nodupes (segmentVP g 0.75 1750 250). Using this g as the first coordinate and the periodic points as second coordinates, plot a small point for each pair (using NetLogo’s plotxy command). Repeat this for enough values of g to ensure an informative bifurcation plot. The following plotBifurcations implementation produces a bifurcation plot very similar to that in the Logistic Growth chapter.

Extending the Logistic Growth Model (NetLogo)

• You may use BehaviorSpace to export the data. (See the Introduction to NetLogo supplement for an introduction to BehaviorSpace.)

Additional Resources

NetLogo’s documentation is very good, but sometimes additional resources can help new users. Recall that the Introduction to NetLogo supplement includes a short introduction to BehaviorSpace. A web search turns up good video discussions of BehaviorSpace. Here are a couple to start with.

• BehaviorSpace in NetLogo by Maureen Psaila-Dombrowski

• Agent-Based Modeling: BehaviorSpace by Bill Rand

• NetLogo’s High Ceiling by Jacob Kelter

The following video discusses NetLogo’s csv extension.

• Netlogo Import and Export Csv Data by Michael Rejtig

LogisticGrowth01 (Python)

Create a new model file named logistic_growth01.py. This file should hold all the code for the LogisticGrowth01 project.

def test_nextVP():
    assert (POP0 == nextVP(0, POP0))
    assert (0 == nextVP(GMAX, 0))
    assert (1 == nextVP(GMAX, 1))

from logistic_growth01 import trajectory
print(*trajectory(), sep='\n')

from logistic_growth01 import trajectory, exportLG
exportLG(trajectory())

from logistic_growth01 import trajectory, plotLG
plotLG(trajectory())

if __name__ == "__main__":
   plotLG(trajectory())

python3 logistic_growth01.py

LogisticGrowth02 (Python)

Begin the LogisticGrowth02 project by refactoring the trajectory computation. As described in The Logistic Growth chapter, implement the new computation as a trajLG02 function. Remember to test your work. A properly implemented trajLG02 function, which takes three arguments, must pass the following test.

def test_trajLG02():
    traj = trajLG02(GMAX, PSS, POP0) #baseline trajectory
    assert (51 == len(traj))
    assert (POP0 == traj[0])
    assert (nextVP(GMAX, POP0/PSS)*PSS == traj[1])
    traj = trajLG02(0.0, PSS, POP0) #zero-growth trajectory
    assert (51 == len(traj))
    assert (POP0 == traj[0])
    assert (POP0 == traj[-1])

Once your trajLG02 passes this test and you are satisfied with your implementation, compare it to the following Python code. This is one possible Python implementation of trajLG02.

With trajLG02 in hand, it is possible to approach the collectTrajectories activity. There are many possible approaches to this activity. In every case, we must produce the three needed trajectories, one for each value of \(g_\text{max}\). Abstractly, this effectively maps a partially applied trajLG02 function over a sequence containing the maximum growth rates. The result of this operation is a sequence of trajectories. Concretely, even though map is a builtin function, a new Python programmer may be more likely to build up a list of trajectories inside a for statement. This is a fine approach, but a more experienced programmer will probably use a list comprehension.

List comprehension is usually the most idiomatic way to map a function across a list of values. The following gmaxExperiment function takes this approach, but it then goes further. Instead of simply returning the list of trajectories, it uses them to produce a Polars dataframe that contains the trajectories as columns. Dataframes are the most common way for data scientists to work with columnar data, and this approach will soon pay dividends.

The gmaxExperiment function takes one argument, a sequence of maximum growth rates. It produces a trajectory for each one of these, using list comprehension. It is possible to produce a dataframe directly from this list of trajectories: Polars will provide each with a default name. Here however, we construct an informative name from each growth rate and provide the dataframe constructor with a list of names by means of its schema argument. We get an immediate payoff from the use of Polars when we attempt to complete the exportTrajectories activity: a Polars dataframe has a write_csv method. All we need to do is provide a filename.

world = World(g=0.03, pop=7.75, pbar=11)
world.runsim()

from matplotlib import pyplot as plt
plt.plot(world.history)
plt.show()