Complex Dynamical Systems

I believe that nothing passes away without leaving a trace, and that every step we take, however small, has significance for our present and our future existence.

—Anton Chekhov, My Life

Not only in research, but also in the everyday world of politics and economics, we would all be better off if more people realized that simple non-linear systems do not necessarily possess simple dynamical properties.

—[May-1976-Nature]

Overview and Objectives

This lecture builds on the programming and data visualization skills introduced in the Exponential Growth lecture. This time, the focal model produces projections for a population whose rate of growth is not constant. The model includes a notion of population saturation, in line with the literature on environmental carrying capacity. The core of the present lecture is the complete development of a simple logistic model of population growth. This is the second and last system-dynamics model in this course.

Before attempting to implement the model of this lecture, skim the entire lecture from start to finish. Be sure to garner a good sense of the conceptual model and of some possible ways to visualize simulation outcomes. Then work through the lecture step by step, implementing the computational model as you go. Continually keep the implementation runnable by temporarily using stubs during model development. Verify each function implementation, either by using provided tests or by writing your own.

Goals and Outcomes

The central goal of this lecture is to enlarge the modeling and simulation toolkit we began developing in the Exponential Growth lecture. Once again, a simple system-dynamics model of population growth provides the backdrop for the introduction of core simulation tools. This time, there is feedback from a stock variable to a flow variable. Specifically, the population growth rate is no longer exogenous: population saturation produces a change in system behavior. A key goal of of this lecture is to demonstrate that, despite the simplicity of the model, this small change can produce radical changes in system behavior. Indeed, for certain parameter values, this simple system engenders complex dynamics.

Prerequisites

Each section of this lecture includes exercises that are structured to foster deeper understanding through active exploration of the core material. These exercises introduce and refine programming skills that future lectures require. The exercises below presume the mastery of a small amount of preliminary material, indicated by the following prerequisites.

State-Dependent Growth Rate

After completing this section, you will be able to do the following.

  • Explain the concept of environmental carrying capacity.

  • Explain the Verhulst-Pearl logistic growth model.

  • Explain the concept of a model parameter.

  • Copy files as a very primitive form of code reuse.

This section introduces the Verhulst-Pearl characterization of population growth. This characterization allows the population growth rate to be state-dependent: the current population growth rate depends on the current level of the population. More precisely, in the present lecture, the current growth rate depends on the relationship between the current population and its sustainable level.

Recall from the Exponential Growth lecture that an evolution rule describes how a dynamical system that is in one state proceeds to its next state. In the exponential-growth population model of the Exponential Growth lecture, the state of the system is the current population. In that model, one year’s population and an exogenous growth rate determine the next year’s population. In terms of the stock and flow concepts of system dynamics, the flow of net population additions changes the stock of people, which is the world population. This lecture endogenizes the population growth rate. Now, the world population (a stock variable) affects the population growth rate (a flow variable).

Conceptually, the growth rate of a population is the difference between the birth rate and the death rate. The conceptual model of the present lecture embeds an implicit model of constraints on population growth. As a population increases, birth rates eventually fall and death rates eventually rise due to environmental constraints. As a result, the growth rate of the population depends on the level of the population. From a modeling perspective, the population growth rate is no longer exogenous. It is now state-dependent.

Constrained Population Growth

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

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 lecture speaks instead of the sustainable population, and later sections link this to the mathematical concept of a steady state.

Mathematical Formulation

The characterization of population change in this lecture is known as the logistic-growth model. To introduce this model, let the annual growth rate of the population depend on the ratio between the current population (\(P_t\)) and the sustainable population (\(P_{ss}\)). The following characterization of logistic growth also introduces a maximum possible growth rate, \(g_\text{max}\). In contrast to the exponential population growth model of the Exponential Growth lecture, this maximum annual growth rate is not the actual annual growth rate. The actual growth rate is not constant; it falls as the population rises.

\begin{equation*} \frac{P_{t+1} - P_{t}}{P_t} = g_\text{max} * (1 - P_t / P_{ss}) \end{equation*}

To simplify notation, define \(p_t\) to be the ratio of the actual population to the sustainable population: \(p_t = P_t /P_{ss}\). The present lecture calls \(p_t\) the relative population. The logistic growth equation may be more simply expressed by using this single variable, as follows.

\begin{align*} \frac{P_{t+1}/P_{ss} - P_{t}/P_{ss}}{P_t/P_{ss}} &= g_\text{max} * (1 - P_t/P_{ss}) \\ \frac{p_{t+1} - p_{t}}{p_t} &= g_\text{max} * (1 - p_t) \end{align*}

The relative population \(p_t\) summarizes the state of this dynamical system at time \(t\). The growth rate of \(p_t\) depends on the level of \(p_t\). That is, the growth rate is state dependent. As environmental constraints become increasingly relevant, growth slows. If they are exceeded, population growth can even become negative.

When the population is small relative to the sustainable population, then \(p_t\) is close to zero. In this setting, the population growth rate approaches its maximum (\(g_\text{max}\)). However, as \(p\) rises, population growth slows. As the population approaches the sustainable population (where \(p=1\)), the population growth rate falls to \(0\).

To produce a simple evolution rule for the logistic-growth model, solve for \(p_{t+1}\) in terms of \(p_t\). This evolution rule provides an algebraic summary of the logistic-growth discrete-time dynamical system, which is the focus of the present lecture.

\begin{equation*} p_{t+1} = p_{t} + g_\text{max} * p_{t} * (1 - p_{t}) \end{equation*}

Logistic-Growth Map

This equation summarizes the relationship between next period’s relative population (\(p_{t+1}\)) and this period’s relative population (\(p_{t}\)). The value of \(p_{t}\) in expression on the right-hand side of this equation determines the value of \(p_{t+1}\) on the left-hand side.

Given any maximum growth rate (\(g_\text{max}\)), a univariate real function can embody the evolution rule for \(p\). This is a discrete-time Verhulst-Pearl growth function, which is more commonly called the logistic-growth function or logistic-growth map. [1] For example, setting \(g_\text{max}\) to \(0.03\) produces the following logistic-growth function. This is a univariate function, because it has a single input. It is a real function because the input and the output are real numbers.

\begin{equation*} p \mapsto p + 0.03 * p * (1 - p) \end{equation*}

The Function Parameter

Recall from the Exponential Growth lecture that the variable \(p\) in the function expression is the function parameter. It is bound by the function expression. A function parameter has meaning only within the function expression, not outside of it. Roughly speaking, function parameters are accessible only locally (i.e., inside the function expression).

For example, the use of a function parameter \(p\) in a different function expression would also be bound by that function expression. The two uses have absolutely no relationship to each other. The use of p in one function expression is completely hidden from and unrelated to any other function expression.

Discrete-Time Dynamics of Logistic Growth

As illustrated in the Exponential Growth lecture, iteration of the evolution operator for a discrete-time dynamical system produces a system trajectory. In this logistic-growth model, the system trajectory is the relative population trajectory.

For example, consider a specific maximum growth rate (\(g_\text{max}=0.03\)) and any initial relative population (\(p_0\)). Starting with \(p_0\), repeated application of the logistic-growth map produces the implied future relative populations. (The corresponding predictions of population levels are \(P_t = P_{ss} p_t\).) Here are the first two iterations along with a general representation of the relationship between relative populations in adjacent years.

\begin{align*} p_{1} &= p_{0} + 0.03 * p_{0} * (1 - p_{0}) \\ p_{2} &= p_{1} + 0.03 * p_{1} * (1 - p_{1}) \\ \vdots \\ p_{t+1} &= p_{t} + 0.03 * p_{t} * (1 - p_{t}) \end{align*}

Computational Logistic-Growth Function

These simple algebraic manipulations underpin equally simple computational manipulations. The first step is to convert the mathematical definition of this logistic-growth map int a computational function. This computational function accepts a single input, which is a real number representing the current relative population. The output is a real number, which represents the next year’s relative population. The following function sketch describes computational function named nextVP that performs this transformation. [2]

function:

nextVP: Real -> Real

function expression:

\(p \mapsto p + 0.03 * p * (1 - p)\)

function parameter:

\(p\), the relative population

summary:

Return the next relative population by incrementing the current relative population (\(p\)) by \(0.03 * p * (1 - p)\).

Create a new project, named LogisticGrowth01. As a crude starting point, you may copy the implementation of the ExponentialGrowth project, since it shares some common features with this new project. (Later lectures explore better approaches to code reuse.) From now on, make changes only to the new project.

Implement the nextVP function, as described in the function sketch. This function works directly with the relative population—the current population relative to the sustainable population. Initially, test nextVP any way you wish. Then, implement tests—perhaps in a testNextVP procedure—and check that nextVP passes your proposed tests.

Model Parameters as Global Constants

The parameters of a simulation model are fixed values that can affect the simulation outcomes. For example, the maximum growth rate is a model parameter for the upcoming logistic growth simulations. In the intial implementation of the nextVP function, this model parameter is simply hard coded. More commonly, the values of model parameters may be set during a simulation model’s set up phase, but they should not be changed by the simulation schedule.

The population growth model introduced in the Exponential Growth lecture was particularly simple. There, population projections depend only on the initial population and the constant population growth rate. That population growth rate in that model is a model parameter; it remains constant during a simulation. It is also the focal parameter: explorations of the model focus on how changes in the population growth rate affect the population projections.

Global constants can be a particularly simple way to handle model parameters. As so often in coding, a choice between approaches involves trade-offs. For initial implementations of a model, this course strongly favors presentational simplicity over computational sophistication. Accordingly, code examples in this course often introduce global constants. This course goes even further in pursuit of simplicity and encourages initial hardcoding of parameters that are not used in experiments. As your programming skills become more sophisticated, you will often wish to explore ways to purge global variables and hard-coded parameters from your code.

For presentational simplicity, the Exponential Growth lecture hardcoded the population growth rate. (See the Explorations for that lecture for an alternative.) This lecture initially continues that practice, but instead of a fixed population growth rate, there now is a fixed maximal population growth rate.

Population projections still depend on the initial population, of course. However, as shown below, population growth now follows a logistic trajectory. The logistic growth model adds a new consideration. Population growth now takes place in an environment characterized by a maximum sustainable population (\(P_{ss}\)). This is another model parameter.

Setting the values of the model parameters is prerequisite to running a simulation. Setting these values is typically part of the setup phase of a simulation model. In this section, the model parameters are P0, Pss, and gmax. Table t:LogisticParameters specifies their values and provides a model specific interpretation. This table also specifies the maximum growth rate as a model parameter, even though the initial model implementation hardcodes that value (in the nextVP function).

Parameters of the Logistic-Growth Model

name

value

meaning

P0

8.1

initial world population (in billions)

Pss

11.5

sustainable population (in billions)

gmax

0.03

maximum annual growth rate

Estimates of the planetary carrying capacity vary substantially. Although \(10\) billion is an often cited estimate, the specification in this lecture is slightly optimistic and raises that by 15% to \(11.5\) billion. [3] This number plays the role the sustainable population (\(P_ss\)) in this lecture. Additionally, assume a maximum population growth rate of 3% per year. As shown below, this produces a plausible initial growth rate.

Global Variables

In addition to setting the values of the model parameters, the setup phase includes the intialization of any global variables. For example, the exponential growth model of the Exponential Growth lecture 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 lecture, 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.

Marking Time

Discrete-time simulations often include a global variable to track the time elapsed. Typically this is just the number of iterations. It is quite common to store this in a global variable named ticks, capturing the idea of the ticking of a simulation clock.

This section correspondingly introduces ticks as a new global variable, which counts the number of simulation steps. It is typically 0 at the end of the setup phase and increments at the end of each iteration of the simulation schedule. This counter does nothing more than keep track of the number of iterations of the simulation schedule. Therefore in this model of population growth, one tick represents the passage of one year. This is the time scale of the model.

Initialization

Recall that simulation models in this course 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 intialize ticks to \(0\) during the setup phase.

Global Variables in the Logistic-Growth Model

name

initial value

meaning

relpop

P0 / Pss

relative population

ticks

0

completed iterations

Modify your LogisticGrowth01 model to declare all the model parameters and global variables and intialize them. This should be completed during the setup phase.

Good practice generally prescribes adding code comments wherever they might prove helpful to a reader. So in your source code, include comments that specify the type and meaning of each model parameter or global variable. Distinguish between Real and Integer types, even if your implementation language does not.

Updating the World Population

Recall from the Exponential Growth lecture that following the setup phase, a typical discrete-time simulation includes an iteration phase. During this phase, the simulation schedule runs repeatedly, until the simulation ends. At the core of this schedule is the model step, which includes the computational evolution rule. (This is the computational implementation of the evolution rule derived from the conceptual model.) Each iteration of the model step changes the model state.

In this section’s logistic-growth population model, the model state is essentially the relative population. The model step uses the nextVP function to update the relative population (code:relpop). Additionally, the model step should increment ticks. This course typically suggests that the core simulation schedule be implemented as a step activity. The following activity sketch is illustrative.

activity:

step

summary:

  • Update relpop to its next value, produced by applying nextVP to the current value.

  • Increment ticks.

Implement the step activity for the logistic-growth population model and add it to the LogisticGrowth01 project.

Producing a Population Trajectory

The two core components of this simulation model are the setup phase and the simulation step. After implementing these in code, it is already possible to run a logistic-growth population simulation. For example, to project relative populations over \(50\) years, simply set up the model and then iterate the simulation step \(50\) times. The following activity sketch describes a runsim01 activity that does this.

activity:

runsim01

summary:

  • For the setup phase, set up the simulation model as above.

  • For each of \(50\) iterations, run the simulation step to produce the next value of the relative population.

That is all well and good, but this runsim activity is not very useful unless supplemented some way to view the population projections. As illustrated in The Exponential Growth lecture, the simplest export of the population data just prints the population values to a computer screen. There are two basic approaches to data output. Either output all the data at once at the end of the simulation, or output it incrementally as the simulation runs. For the simulations in this course, the two approaches are generally equally satisfactory.

Implement the runsim01 activity, described above. Then augment it to print the population trajectory.

Improved Data Handling

Although it may be possible to copy printed output and paste it into a spreadsheet, printing the population projections to a computer screen is neither elegant nor convenient. This section therefore explores a different, more useful approach: exporting the data to an output file. This requires a rather small change in the simulation. Instead of printing population projections to a computer screen, running the simulation writes them to a file.

The Output File

Export of the simulation data to a file involves two conceptual steps: setting up an output file, and writing simulation data to this output file. These steps differ substantially across languages and toolkits, and they also differ between incremental and end-of-simulation bulk export strategies. Conveniently, simulation toolkits often offer facilities for data export that require very little effort from the user. In contrast, some programming languages demand a substantial understanding of file handling.

Create an output folder named out to hold the simulation output. (This lecture assumes that that the out folder is immediately below the folder containing the LogisticGrowth01 source code.) Augment the LogisticGrowth01 project as needed to export the population projections to a file name LogisticGrowth01.csv in the out folder. For simplicity, just hardcode the file path.

File Handling for Data Export

The following more detailed discussion of data export may be skipped for the moment and returned to as needed. This discussion addresses two strategies for data export: incremental export, and end-of-simulation bulk export. Planning for either requires choosing a file format. A particularly simple plain text format writes just one number per line. Spreadsheet and statistical applications can readily import this format.

Bulk Export

In order to perform an end-of-simulation bulk export, the simulation data much be accumulated by the simulation program. In the present case, there is a single value that is tracked: the population. Collecting this data as the simulation runs normally involves appending a new population value each period to some kind of mutable sequence. After the last iteration, the collected data includes all of the needed projections. Before exiting, the program should write this data to file. This becomes a wrap-up phase for the simulation.

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.

This file handling exercise is optional. You may use BehaviorSpace to export your data.

Create a setupFile activity, as described above. Create a recordStep activity, as described above. For simplicity, simply hardcode the file path (e.g., "out/LogisticGrowth01.csv").

Running the Simulation

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 initializes the global variables. 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. The following activity sketch describes such an activity, and the subsequent exercise elaborates on its use.

activity:

runsim

summary:

  • Execute the setup phase of the simulation model.

  • Repeatedly (\(50\) times) execute the iteration phase of the simulation model.

Use the model parameters listed in Table t:LogisticParameters. Implement the runsim activity described above. Produce 50-years of population projections, exported to a file. Load the exported simulation data into a spreadsheet and produce a run-sequence line plot of the data.

Comparing Population-Projection Plots

Plots often make it easy to see patterns in data that might be challenging to detect in raw numbers. Figure f:CompareExpLog presents a chart based on the data generated by this new simulation along with the simulation from the Exponential Growth lecture. The difference between the projection based on exponential growth and the projection based on logistic growth is immediately evident in the plot. 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.

Compare the exponential and logistic growth population projections.

Population Projections (Exponential vs Logistic Growth)

Initial growth rates begin at about 1% annually. This is the constant growth rate in the exponential-growth model. The logistic-growth model set a maximum growth rate of 3% annually, but the actual growth rate falls as the population approaches its sustainable level.

A Simulation Experiment for the Logistic-Growth Model

After completing this section, you will be able to do the following.

  • Compare the mathematical concept of a free varaible to the computational concept of a nonlocal variable.

  • Explain the concept of a simulation experiment.

  • Experiment with a simulation by varying a model parameter.

  • Create useful data visualizations to support a simulation experiment.

The primary goal of this section is to explore the idea of a simulation experiment. A secondary goal is to elucidate some implications of logistic growth. The project supporting these goals is a very basic experiment with the logistic-growth population model. (Later lectures consider more sophisticated experiments.)

Role of Maximum Growth Rate

Running the logistic-growth simulation model of the previous section produces a set of population projections. These eventually diverge substantially from the exponential growth projections of the Exponential Growth lecture. The 50-year forecast from the logistic-growth population model is much lower than the exponential-growth forecast. The new forecast appears more realistic. [4]

It is natural to wonder how sensitive these population projections are to a key model parameter: the maximum population growth rate (\(g_\text{max}\)). This constitutes a research question that may be addressed by means of a simulation experiment. This section provides a first approach to the question of how the population projections change in response to changes in the maximum-growth-rate paramater. The method is to produce the population trajectories for several different maximum growth rates and then compare them.

Primitive Experimentation

As a first approach to experimentation in this model, try out the following error prone and tedious strategy. Repeatedly change the definition of nextVP, rerun the simulation, and recreate the associated spreadsheet plot. The point of this exercise is to gain insight into the limits of such an approach to simulation experiments.

Repeatedly change the nextVP function so that the hardcoded maximum annual growth rate is successively \(1\%\), \(3\%\), and \(6\%\). For each value, rerun the simulation and use a spreadsheet to illustrate the resulting population projections. Remember to reset the maximum growth rate to \(3\%\) when you are done with this exercise.

Some Problems with Primitive Experimentation

The maximum growth rate is the focal parameter for this experiment: we wish to change it and discover how this changes the simulation outcomes. Unfortunately, it is hardcoded in the source code. Repeatedly changing a hardcoded parameter by hand is a tedious approach to experimentation. It requires revisiting the code, finding the hardcoded value, and altering it. This practice becomes more tedious with more complicated models or with more focal parameters. It is also error prone. It can be hard to remember the last value set, and there is no hint in the output file of the parameterization used to produce it.

Despite such drawbacks, such primitive experimentation may occasionally prove useful very early in the modeling process. Even so, it is important to move quickly beyond it. This course introduces several better approaches, and the present lecture focuses on a small first improvement: eliminating the repeated changes to a hardcoded value that is buried inside a function definition. To accomplish this, the model code will introduce a name for the maximum growth rate rather than just working with its value.

Hardcoding vs Explicit Parameterization

A model parameterization assigns values to all the model parameters. During a simulation run, the model parameters retain the same fixed values. In a simple system-dynamics model without exogenous inputs, a complete model parameterization fully determines the behavior of the simulation. A simulation experiment varies the model parameters across simulation runs in order to explore the resulting changes in system behavior.

In the Exponential Growth lecture, the definition of the exponential growth function included a hardcoded 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). So far, this lecture has adopted the same approach: the hardcoded maximum growth rate ensured that nextVP was a pure function.

For the current simulation experiment, however, hardcoding 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 setting it to varied values. This is particularly true when it is buried in a function definition, somewhere in the source code. It becomes even more problematic if a hardcoded parameter appears in multiple locations in the program code, so that changing it requires changing the code in multiple locations. It is wise to assume that any attempt to coordinate such code changes eventually will fail.

On the one hand, pure functions are a good thing: they are easy to understand and have completely predictable behavior. On the other hand, it is a good idea for our function definitions to avoid hardcoding values that might frequently change. This lecture explores one possible resolution of this tension: we accept that functions definitions sometimes may usefully refer to nonlocal variables. Later lectures explore other approaches.

Nonlocal Variables

When programmers say that a function definition depends on a nonlocal variable, they mean that there is an outside influence on the values returned by the function. The value of a nonlocal variable exists outside of the function definition. It is part of the environment in which the function is executed.

The definition of a pure function provides a complete characterization of the relationship between its inputs and its outputs. In contrast, when a function depends on a nonlocal variable, the relationship between a function’s explicit inputs and its output is unknow 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 a nonlocal variable in the conduct of simulation experiments.

(Advanced) Free Variables

Nonlocal variables in a computational function resemble free variables in a mathematical function. For example, the variable \(g_\text{max}\) is free in the following mathematical function. It represents a value that is not supplied as a function argument but is still needed in order to specify the function’s numerical value. This free variable has meaning outside the function expression. In contrast, the variable \(p\) is bound in the following mathematical function. It has no meaning outside the function expression.

\begin{equation*} p \mapsto p + g_\text{max} * p * (1 - p) \end{equation*}

A mathematician would say that this function definition includes \(g_\text{max}\) as a free variable, because \(g_\text{max}\) is not bound by the function definition. This course freely uses this mathematical terminology. A computer scientist is more likely to say that an implementation of this function is not pure unless it replaces \(g_\text{max}\) with a number. Otherwise, gmax is nonlocal to the function definition.

The important point is that, based on this function definition alone, the value of \(g_\text{max}\) is unknown. This implies that the output of the function is unknown, even when we precisely know the input. Roughly speaking, in order for the output of a function to be completely determined by its arguments, the function definition must contain no free variables. That is, the definition must be completely closed to outside influences.

For example, to be closed in this way, the definition of a numerical function must completely specify the numerical relationship between the inputs and the outputs. The original logistic-growth function definition, \(p \mapsto p + 0.03 * p * (1 - p)\), is closed; it contains no free variables. In contrast, the function definition \(p \mapsto p + g_\text{max} * p * (1 - p)\) is open to outside influences. It contains the free variable \(g_\text{max}\).

Unlike the function parameter \(p\), \(g_\text{max}\) is not bound in this function definition. This is evident because the variable \(g_\text{max}\) does not occur to the left of the mapsto arrow; it is not a function parameter. Unlike the function parameter \(p\), the variable \(g_\text{max}\) is not bound in this function definition. This means that \(g_\text{max}\) is a free variable in this function definition. A function definition containing a free variable is called open, because the return value is open to influences from outside of the function definition.

Roughly speaking, in order for the output of a function to be completely determined by its arguments, the function definition must contain no free variables. That is, the definition must be completely closed to outside influences. For example, to be closed in this way, the definition of a numerical function must completely specify the numerical relationship between the inputs and the outputs.

However, an advantage of this open logistic-growth function is that it provides a general representation of the functional form for logistic growth. Many particular univariate functions have this functional form—one for each particular maximum growth rate. Working with free variables can therefore be useful. For example, working with the open definition of logistic growth may enable us to reach some conclusions about population growth that will be true for any value of the free variable math:g_text{max}.

A Revised nextVP Function

In order to experiment with changes in the maximum growth rate, this section changes the computational model to explicitly parameterize the maximum growth rate (gmax). This parameter will be a nonlocal influence on the values produced by the nextVP function. Changes in the value of gmax therefore affect our population projections. The simulation experiment of this section explores these effects.

The model changes required to parameterize the maximum growth rate are small, but they fundamentally change the logistic-growth function. The original version of nextVP hardcoded the maximum growth rate as \(0.03\), which produced a pure function. The new version of nextVP replaces that hardcoded \(0.03\) with the model parameter gmax. This means the behavior of the function remains unknown until we specify the value of the model parameter gmax.

When computing numerical results, the nextVP function now depends on the value assigned to gmax. Each value assigned to gmax effectively produces a different specific logistic-growth function, which affects the entire trajectory of the logistic-growth population simulation. Since the value of gmax cannot be determined from the function definition alone, gmax is nonlocal in the function definition. (If gmax is visible to the entire model, its scope is not just nonlocal, it is global.) The behavior of the function now depends not only on the explicit function definition but also on the context in which it runs.

Create a new project, named LogisticGrowth02. (You may begin by making a copy of the LogisticGrowth01 project.) Introduce a global constant, gmax. For now, the setup phase of your simulation should ensure that gmax has 0.03. Replace the definition of nextVP with a new version, which depends on this model parameter. (The maximum growth rate is no longer hardcoded in nextVP.)

Context

The new nextVP function cannot produce numerical results outputs except in a context that specifies a value for \(g_\text{max}\). This new function is not closed to outside influences; it depends on the model parameter \(g_\text{max}\). Many specific functions have the logistic-growth functional form—one for each particular \(g_\text{max}\). Parameterizing the maximum growth rate produces an abstract representation of all of these particular functions.

Each simulation run with a model depends on the values given to the model parameters, which can be different for every simulation run. The new function definition allows us to easily experiment with changes in the maximum growth rate. What is more, this approach required minimal tinkering with the original univariate function definition. This is a particularly simple approach to supporting a simulation experiment.

Experimenting with the Maximum Growth Rate

Parameterization of the maximum growth rate makes the increasing complexity of the model a little bit more manageable: there is a clearly defined unique location for any code changes. A change in the maximum growth rate produces a change in the behavior of nextVP. Take advantage of this to do a simple experiment.

Run the LogisticGrowth02 simulation for three values of gmax: 1%, 3%, and 6%. Read the resulting population projections into a spreadsheet and produce a chart comparing the three trajectories.

If you use BehaviorSpace, first move gmax to a slider. The create a new experiment named lg02, and replace the Vary variables contents with ["gmax" 0.01 0.03 0.06]. Recall that the BehaviorSpace Discussion on Canvas includes this useful video.

Renaming 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 course 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.

In the LogisticGrowth02 project, create an outfile variable that holds the path to the output file. In the setup phase (e.g., in the setupFile activity), set the value of outfile to LogisticGrowth-gmax_{gmaxpct}.csv, where {gmaxpct} represents gmax as a percentage. (Multiply gmax by 100.) For example, when gmax is 0.03, the filename should be LogisticGrowth-gmax_3.csv. Modify the data export code to use this new filename.

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.

Compare logistic growth projections for various maximum growth rates.

Logistic Population Projections (Various \(g_\text{max}\))

Short-Run vs Long-Run Population Projections

In this sense sense, population projections are sensitive to these parameter changes. These changes in the maximum growth rate have fairly large effects on the population projections. In another sense, however, the simulation outcomes are eventually invariant to these changes in the parameterization: the population always converges to the sustainable population. Changing the maximum growth changes the speed of convergence but not the sustainable population. This implies that, for the values of gmax considered, the long-run population projection does not depend on the maximum growth rate.

GUI Experimentation

  • Create a slider widget for interactive experimentation.

  • Create a button widget for simulation control.

  • Create an animated time-series plot.

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.

Widgets for Simulation Control

So far, this project has utilized a command line interface in order to run the population simulations. This is quite flexible, but for interactive experimentation it is not fully satisfactory. For one thing, it may involve tedious and error prone typing at the command line. Additionally, it complicates sharing of the model: command-line usage may appear obscure to new users of the model. A standard remedy for these problems is to add a GUI interface to the model.

Clickable buttons are a particularly common GUI widget. (A GUI widget is a graphically presented control or display element, such as a control button, a slider, a text box, or a plot element.) Models built with simulation toolkits often provide buttons that allow even a new user to easily set up and run the simulation. Correspondingly, add a button for simulation control to the LogisticGrowth02 project.

Add a GUI button labeled Run Sim to your LogisticGrowth02 project. For the moment, this Run Sim button should just execute the runsim01 activity. Once you have created your button, use it to run your simulation. This should print the same propulation projections as before.

Sliders

Watch this video on slider creation.

One popular approach to interactive parameter manipulation, particularly in the early phases of simulation modeling, is the addition of control widgets. A slider is a control widget that provides a simple graphical interface for setting a variable’s value. It is a natural choice when we wish to experiment with a parameter that can take on a range of numerical values. A slider typically specifies a default value for a variable and then allows a user to reset that value within a constrained range.

The next goal of this section is to demonstrate the use of a GUI slider for interactive model experimentation. Accordingly, the next exercise adds to the model GUI a slider that controls the maximum growth rate. Remember that model parameters should not change during the iteration phase of the simulation. So, use this slider only to influence the model setup, which must take place after the slider is set.

Add a GUI slider to your LogisticGrowth02 project. This should control the value of the maximum growth rate with the new global variable (gmax). Design the slider so that it can take on at least the values of 1%, 3% (the default), and 6%. (A later section considers additional values.)

Time-Series Plot

The Exponential Growth lecture illustrated the use of dynamic plots for the display of simulation results. This section continues that practice and then discusses widgets for GUI experimentation with a simulation model.

Prepare for this by adding a time-series plot of the population projections to your LogisticGrowth02 model. Creation of the time-series plot should closely follow the discussion in the Exponential Growth lecture. Plotting must incorporate one substantive difference from that previous work: since the core simulation produces a trajectory for the relative population, you must transform the data to plot the population levels. Rely on the previous discussion of the runsim01 activity to guide this transformation.

First, review the discussion of plotting in the the Exponential Growth lecture. Then, add a time-series plot to the LogisticGrowth02 project. This plot should display the annual population projections produced by the population simulation. (If possible, this should be a dynamic plot.) The title of this plot should be "Time Series".

To populate the plot, implement a runsimG activity that modifies your runsim activity by adding plotting. Remember to include plotting of the initial model state in the setup phase. Use this runsimG activity to plot the annual population projections for \(50\) years.

Refinements

In anticipation of future models, add a few refinements to the existing LogisticGrowth02 code. Isolate the plot setup, and bundle the new simulation schedule.

Implement an activity that sets up the plot; call it setupPlotTS. This activity should plot the initial population and, if possible, set the axes limits and labels.

As a minor convenience, package the entire simulation schedule (including plotting) into a new go activity. (This will execute your step activity.) Refactor your runsimG activity to repeatedly execute this new go activity.

GUI Exploration of the Logistic Growth Model

During the early stages of simulation modeling, interactive experimentation is common. The GUI of the LogisticGrowth02 project is now rich enough to support this. Set the new slider for gmax to the default values of \(0.03\), and use the Run Sim button to run the simulation. Next, use the slider for gmax to double the maximum growth rate to \(0.06\). The behavior of the population is very similar, but convergence to the \(P_{ss}\) is faster. Also use the slider for gmax to cut the maximum growth rate to \(0.01\). The behavior of the population is again very similar, but convergence to the \(P_{ss}\) is slower. The results should not be surprising: they resemble Figure f:CompareExpLog. This rapid process of interactive experimentation illustrates the usefulness of plots, buttons, and sliders for model exploration.

(Advanced) Complex Logistic Growth

This section includes some more advanced material that is not part of the assignment. Readers who wish to implement the examples in this section will need to know how to use their simulation toolkit to plot a function.

This section further explores the logistic-growth function by considering radically different values of gmax. These values do not remotely describe human population growth, which correspondingly is not the focus of the present section. Rather, this section explores a few challenging mathematical concepts and more advanced programming skills. Some readers may prefer to skip this section for now and return after mastering additional models in this course.

Revised Plot Code

Revised Slider

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

Create a new LogisticGrowth03 model by copying the LogisticGrowth02 model. Do not make any further changes to the old model. Alter the slider in the new model to extend from 1% to 300%.

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.

Change the commands run by the Set Up and Go 50 buttons. Setting up the simulation now includes setting up a population plot. Since the plot setup includes plotting the initial state of the model, it should come after initialization of the model. The Go 50 button should not only run the simulation step but should also update the time-series plot.

NetLogo builds in some nice refinements. Examine the button creation dialog, and note the checkbox labeled “Disable until ticks start”. Since it does not make sense to attempt a simulation step until the model setup is complete, check this box for the Go 50 button.

Visualizing Logistic Population Growth

Steady State

Begin with the concept of a steady state. A steady state of a dynamical system is a state that, in accordance with the evolution rule of the system, does not change over time. If a dynamical system is in a steady state, it stays there.

In the previous section, the sustainable population determined the long-run outcomes of the logistic-growth simulations. Furthermore, \(P_{ss}\) is a steady-state population. That is, a population equal to the \(P_{ss}\) will not change.

Steady States as Fixed-Points

From a dynamical point of view, it is natural to speak of steady states of a dynamical system. A related concept from mathematics is a fixed point of a function. When application of a function to a value reproduces that value, we say the value is a fixed point of the function. So a steady state is a fixed point of the evolution rule.

Consider the logistic-growth function: \(p \mapsto p + g_\text{max} p (1 - p)\). At a steady state (i.e., unchanging) population, the output value of this function is the same as the input value. We may accordingly ask, for what values of \(p\) is this the case?

Simple algebra can find the possible steady states of the logistic-growth model. The are the solutions of the following equation in \(p\).

\begin{equation*} p = p + g_\text{max} * p * (1 - p) \end{equation*}

This just states the requirement that the function produce a value equal to its input. By substitution, we see that both \(0\) and \(1\) are fixed points of the logistic-growth function. This means that they are steady states of our dynamical system. Recall that a relative pouplation of \(1\) corresponds to a population that equals the sustainable population (\(P_{ss}\)). As we might guess intuitively, a population of \(0\) is also a steady state.

These two fixed points have very different properties. The population simulations of the previous section always converge to a relative population of \(1\); never to 0. The value \(1\) is a stable fixed point of our evolution rule: when the population is near the sustainable population, it approaches it. Equivalently, we say that the value 1 is an attracting fixed point, or attractor, of this dynamical system. 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 this dynamical system.

Picturing Fixed Points

A common visualization of the fixed points of a univariate real function considers a function plot in a two-dimensional coordinate system (i.e., the Cartesian plane). After adding a plot of the \(45^\circ\) line, we may visually determine points where the function plot crosses this line. These points of intersection identify the fixed points of the function.

Create a "Logistic Function" chart for the LogisticGrowth01 model. Create a setupFuncPlot activity that plots a \(45^\circ\) line and the nextVP function in this chart. Consider a plot domain of \([0.0\,..\,1.5]\), but plot points only when the function value is positive. Explore the implications of the parameter gmax by repeatedly changing its value and examining the resulting function plot. Once again consider values of \(0.03\), \(0.06\), and \(0.12\) for gmax. Are the fixed points sensitive to the value of gmax?

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\).

Fixed points of logistic-growth function.

Fixed Points of nextVP (\(g=1.6\))

New Short-Run Behavior at Higher Growth Rates

Slider for the Logistic-Growth Population Model

name

default value

range

increment

gmax

\(0.03\)

\([0\;..\;3]\)

0.01

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.

As a small convenience for the upcoming experiments, make two changes to the model GUI. In addition to setting up the simulation, the Set Up button should now set up both charts. Furthermore, add a Go 5 button to the model, which runs go a total of \(5\) times. Setup the model and run the simulation for \(5\) iterations.

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).

Time Series from logistic growth function iteration when :code:`gmax` is :math:`1.6`.

Time Series from nextVP Iteration (\(g=1.6\))

New Long-Run Behavior at Higher gmax

Further increases in gmax produce new surprises in system behavior. To illustrate, use the slider to set gmax to \(2.2\) and run the simulation. Once again the simulated sequence is not monotone. This time however, the fluctuations do not appear to die out. In fact, as illustrated by Figure f:vpts220, increasing the number of iterations makes it clear the system trajectory quickly settles on two values and alternates between them.

Time Series from logistic growth function iteration when :code:`gmax` is :math:`2.2`.

Time Series from nextVP Iteration (\(g=2.2\))

Periodic Attractors

This is an example of a periodic attractor, and the two values are called periodic points. More specifically, this is a period-2 cycle. This is an example of a limit cycle: the cycle is the limiting behavior of the trajectory.

Lost Steady State?

When gmax is \(2.2\), the trajectory displays a limit cycle. Limit cycles are a radical change in system behavior: the trajectory shows no tendency to approach the steady state over time. In contrast, when gmax is \(1.6\), the trajectory quickly approaches a stable steady-state (Pss). Yet in both cases, the steady-state algebra suggests that the steady states of the system are unchanged. How can we explain these very different trajectories?

Examination of the function plot for logistic-growth when gmax is \(2.2\) confirms there is still have a fixed point at \(1\). 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. [5] The sustainable population is still a steady state for the system, but it is no longer a stable steady state.

Periodic Attractors as Fixed Points

So the steep slope explains the lack of convergence, but what explains the periodic behavior? Consider the periodic points that are produced when gmax is \(2.2\). Suppose one starts with one of these points, applies the logistic-growth function to it, and then applies the logistic-growth function to the result. This reproduces the original point. That is, this point is a fixed point of the twice iterated function.

As a convenient notation, for any function \(f\), represent the twice iterated function as \(f^{\circ 2}\). Again using the graphical approach to discovering fixed points, plot \(f^{\circ 2}\) and search for intersections with the \(45^\circ\) line. Since any fixed point of \(f\) must be a fixed point of \(f^{\circ 2}\), the original steady-state points are rediscovered. For example, when \(f = p \mapsto p + 2.2 p (1 - p)\), we find \(f^{\circ 2}\) still has fixed points at \(0\) and \(1\). But in addition, there are now \(2\) more intersections: the two periodic points (\(0.75\) and \(1.16\)). Figure f:log02g220 illustrates these four fixed points.

The four fixed points of the twice-iterated logistic-growth function when the maximum growth rate is 2.2.

Four Fixed Points of Twice-Iterated Logistic Growth (\(g=2.20\))

Longer Cycles

Our discovery of period-2 cycles shows that our iteration of our simple logistic-growth function can generate surpisingly complicated dynamics. However, there is more to come. To illustrate, set \(g=2.5\) and rerun the simulation. Start with \(5\) iterations, and increase the number of iterations until you see a pattern. As illustrated by Figure f:vpts250, the resulting plots are a bit harder to scrutinize. However, it seems that the time-series plot may be oscillating across \(4\) different values. This is in fact the case.

Time Series from logistic growth function iteration when :code:`gmax` is :math:`2.5`.

Time Series from nextVP Iteration (\(g=2.5\))

Period Doubling

Recall that when gmax is \(1.6\), the simulation trajectory converges to a single value (\(1.0\)). In the long-run, only this value is relevant. When gmax is \(2.2\), the simulation trajectory does not converge to a single value (\(1.0\)). Instead the trajectory is attracted to a period-2 limit cycle. Two points have long-run relevance for this system, which alternates between them. When gmax is \(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 will continue to result from increases in gmax.

The answer is yes, and with faster and faster period doubling. (This is called a period-doubling cascade.) Unfortunately, it becomes harder and harder to detect these oscillations visually. For example, seting \(g=2.55\) produces a period-8 limit cycle, but as seen in Figure f:vpts255, some of the points are too close together to distinguish visually.

Time Series from logistic growth function iteration when :code:`gmax` is :math:`2.55`.

Time Series from nextVP Iteration (\(g=2.55\))

Burn-In Iterations

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. As a very simple example, even gmax is 1.6, the long-run behavior of the trajectory is not immediately evident. To deal with this problem, we often discard initial iterations of a dynamic system in order to get a sense of the behavior it ultimately settles into. The number of discarded iterations is called the burn-in period.

With infinitely precise computations, the actual trajectory would approach the limit cycle over time. Unless the initial point is a periodic point, this approach continues forever. However, computers have computational limitations, so an actual cycle between two numbers emerges after a finite number of iterations. This suggests a strategy for numerically approximating periodic points: run the simulation for a while, and examine only the values towards the end of the trajectory. discard the first part of the trajectory, and then look at the values.

Construct a trajectory of length \(500\) for each of the folloing values of gmax: 1.60, 2.20, 2.50, and 2.55. You may create a function trajectoryVP to produce the trajectory. In each case determine the number of unique values in the final \(100\) iterations.

Bifurcation Diagram

A bifuraction is the doubling of the period of the limit cycle. A bifurcation diagram is popular way to illustrate the effects of changing the parameter gmax. For each value of gmax, it plots the limit 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.

Brute Force Computation of Limit Points

A typical algorithm for brute-force computation of limit points is as follows.

function:

limitPoints: (Real, Real, Integer, Integer) -> List[Real]

function parameters:

  • gmax, the maximum growth rate

  • p0, the initial value of \(p\)

  • n, the number of recorded iterations

  • b, the number of burn-in iterations

summary:

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

  • Let \(p_s = f^{\circ b}(p_0)\), discarding the burn-in iterations.

  • Collect the results of the next \(n\) function iterations.

  • Discard any duplicate results.

  • Return the unique results (i.e., the periodic points).

There are a number of problems with a brute force approach. 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 number of subsequent iterations? Nevertheless, it is by far the most common approach to create a bifurcation diagram.

Brute Force Computation of a Bifurcation Diagram

With a limitPoints function in hand, we can choose a collection of growth rates to examine, and for each growth rate plot all its limit points. A typical result is the following figure.

Limit points of iterated logistic-growth function for various growth rates.

Periodic Points (Iterated Logistic Growth)

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\), we encounter many situations 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 many implications that have spawned a very large literature. Of particular interest, 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 another way, we might say that long-run prediction becomes infeasible.

Summary and Conclusions

Conclusion

Throughout this course we will be repeating these 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 lecture illustrated these components of simulation modeling by developing a simple model of logistic population growth.

Logistic Growth vs Exponential Growth

The logistic-growth model incorporates the concept of a sustainable population. This may be linked to the environmental carrying capacity, which is the population that the the available habitat can sustain indefinitely. This simple increase in realism leads to very different population projections than the exponential-growth model.

Implications for Social-Science Modeling

The discovery of chaotic dynamics in even extremely simple dynamical systems—such as the logistic growth model—can be seen as a negative result for the feasiblity of prediction in the social sciences. Social interactions are inherently nonlinear and dynamic. They are also much more complicated than our simple model of logistic population growth.

Exploration, Resources, and References

Logistic-Growth Explorations

  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. 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.)

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

  4. 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.

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

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

  7. Instead of hardcoding the output file name, introduce it as a global constant (say, outputFile). Change setupFile and recordStep accordingly.

  8. Modify setupFile to write a header line, say population, on the first line of the output file.

  9. (Advanced) Refactor recordStep to take one argument, the current population. That is, instead of having this activity depend on global variables to compute the population, provide a computed value each time it is called.

  10. (Advanced) Instead of hardcoding the number of iterations, make it a model parameter (say, nIter). Control its value with a slider in the model’s GUI.

  11. (Advanced) After comparing to the logistic-growth map in the present lecture, redefine the exponential-growth function in the Exponential Growth lecture to similarly be bivariate.

    Redefine nextVP to have two explicit inputs, gmax and p. This implements the bivariate real function \(\{ g,p \} \mapsto p + g p (1 - p)\). As part of the simulation setup, create a population transition function by partially applying nextVP to the maximum growth rate parameter. Use this new function in your simulation schedule.

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

  13. (Advanced) When parameterizing the name of the output with the growth rate, pad the growth rate with zeros so that the resulting stringified number requires three total characters. For example, if gmax is 0.03, then set the filename to LogisticGrowth-gmax_003.csv, but if gmax is 1.20, then set the filename to LogisticGrowth-gmax_120.csv.

  14. (Advance) The core model description in this lecture includes two global variables. Refactor the model so that it does not need these global variables.

Resources

Population Models and History

The population models in this course use a system dynamics approach rather than an agent-based approach. [Geard_etal-2013-JASSS] offer an agent-based analysis, with an emphasis on demographic structure.

The American Museum of natural history describes the history of human populations in a video entitled Human Population Through Time. For an introductory discussion of the concept of carrying capacity, see the Wikipedia entry. [Sayre-2008-AAAG] discusses varied uses of the term and their historical evolution, and [Hui-2006-EcolMod] clarifies the distinction between population equilibrium and carrying capacity. [Arrow.etal-1995-EcoEcon] and [Niccolucci.etal-2007-EcolEcon] relate carrying capacity and environmental sustainability. [Wilson-2003-Knopf] offers a readable discussion of the sustainability of human populations.

The Logistic Map

Wikipedia includes helpful introductions to the Logistic Map and step plots. [Feldman-2012-OxfordUP] provides a more detailed yet extremely accesible introduction to both topics, and his chapter 7 discusses population models and a logistic function. [May-1974-Science] and [May-1976-Nature] provide early discussions of chaotic dynamics in very simple nonlinear systems.

Chapter 3 of [Kingsland-1995-UChicago] provides a detailed history of the introduction of the logistic hypothesis into research on population growth. [Rosling.Rosling.Rosling-2018-Flatiron] provide a data-rich discussion of the importance of such considerations to human population forecasts. The United Nations Department of Economic and Social Affairs provides world population data and projections, along with helpful charts displaying some of this data. [Cramer-2002-WP] explores the historical development of logistic regression, with connections to the early work of [Verhulst-1838-CMP].

Cobweb Charts (Step Plots)

A popular alternative visualization of these periodic outcomes is provided by a cobweb chart, which includes a type of step plot. Given a sequence \((x_0,x_1,x_2,\ldots,x_n)\), a cobweb chart plots stepped line segments between the points \((x_0,x_1)\), \((x_1,x_1)\), \((x_1,x_2)\), and so on. Conventionally, cobweb chart plots additionally start and end on the \(45^\circ\) line by adding the point \((x_0,x_0)\) at the start and \((x_n,x_n)\) at the end.

Draw two segments each iteration. Given the value \(p_0\) we can produce the value \(p_1=f[p]\), which allows us to draw segments from \((p_0,p_0)\) to \((p_0,p_1)\) and then from \((p_0,p_1)\) to \((p_1,p_1)\). In order to do this however, we need access to both values. In the code, use two global variables (named p and plag) to store these values. Change the project code to include p and now plag as global variables. Each iteration, update both of these globals before you update your plot.

To support this exercise, introduce a new global variable named lagrelpop, which holds the lagged value of relpop. The step activity should now set lagrelpop to relpop before updating relpop. Overlay a step plot on your function plot, producing a cobweb chart. (You may add plotting an initial point to the plot setup.) After each step, the plot update should plot the new point of the function graph and then a new point of the \(45^\circ\) line.

Once the step-plot code is functional, re-examine our last experiment. Do 5 iterations at a time and watch how the plot changes. After about \(20\) function iterations, there are now more visible changes in the step plot.

To pursue this observation a bit further, without clearing the global variables, clear the plots and then plot a few more iterations. In the time-series plot, this renders the period-2 cycle readily visible as an oscillation between two values. In the step plot, we see an unchanging square, reflecting this oscillation between the two period-2 points.

References

[Arrow.etal-1995-EcoEcon]

Arrow, Kenneth, et al. (1995) Economic Growth, Carrying Capacity, and the Environment. Ecological Economics 15, 91--95. http://www.sciencedirect.com/science/article/pii/0921800995000593

[Cramer-2002-WP]

Cramer, J.S. (2002) "The Origins of Logistic Regression". Tinbergen Institute Working Paper 2002-119/4. https://papers.tinbergen.nl/02119.pdf,http://dx.doi.org/10.2139/ssrn.360300

[Feldman-2012-OxfordUP]

Feldman, David P. (2012) Chaos and Fractals: An Elementary Introduction. : Oxford University Press.

[Geard_etal-2013-JASSS]

Geard, Nicholas, et al. (2013) Synthetic Population Dynamics: A Model of Household Demography. Journal of Artificial Societies and Social Simulation 16, 8. http://jasss.soc.surrey.ac.uk/16/1/8.html

[Hui-2006-EcolMod]

Hui, Cang. (2006) Carrying Capacity, Population Equilibrium, and Environment's Maximal Load. Ecological Modelling 192, 317--320. http://www.sciencedirect.com/science/article/pii/S0304380005003339

[Kingsland-1995-UChicago]

Kingsland, S.E. (1995) Modeling Nature. : University of Chicago Press. https://press.uchicago.edu/ucp/books/book/chicago/M/bo3630803.html

[May-1974-Science]

May, Robert M. (1974) Biological Populations with Nonoverlapping Generations: Stable Points, Stable Cycles, and Chaos. Science 186, 645--647. https://science.sciencemag.org/content/186/4164/645

[May-1976-Nature]

May, Robert M. (1976) Simple Mathematical Models with Very Complicated Dynamics. Nature 261, 459--67.

[Niccolucci.etal-2007-EcolEcon]

Niccolucci, Valentina, Federico M. Pulselli, and Enzo Tiezzi. (2007) Strengthening the Threshold Hypothesis: Economic and Biophysical Limits to Growth. Ecological Economics 60, 667--672. http://www.sciencedirect.com/science/article/pii/S0921800906005416

[Rosling.Rosling.Rosling-2018-Flatiron]

Rosling, Hans, Anna Rosling Roennlund, and Ola Rosling. (2018) Factfulness: Ten Reasons We're Wrong About the World---and Why Things Are Better Than You Think. : Flatiron Books.

[Sayre-2008-AAAG]

Sayre, Nathan F. (2008) The Genesis, History, and Limits of Carrying Capacity. Annals of the Association of American Geographers 98, 120-134. https://doi.org/10.1080/00045600701734356

[Verhulst-1838-CMP]

Verhulst, P. (1838) Notice sur la Loi que la Population Poursuit dans son Accroissement. Correspondance Mathematique et Physique 10, 113--121.

[Verhulst-1845-NMAR]

Verhulst, Pierre-Fran,cois. (1845) Recherches Mathematiques sur la Loi d'Accroissement de la Population. Nouveaux Memoires de l'Academie Royale des Sciences et Belles-Lettres de Bruxelles 18, 1--41. http://gdz.sub.uni-goettingen.de/dms/load/img/?PPN=PPN129323640_0018&DMDID=dmdlog7

[Wilson-2003-Knopf]

Wilson, Edward O. (2003) The Future of Life. : Knopf Doubleday Publishing Group.

Appendix: Logistic Growth (Details and Hints)

Logistic Population Growth in Python

Getting Rid of Globals

Python models often avoid introducing the global variables that are natural in NetLogo. One common approach to this is to introduce a World class that is able to initialize and run the simulation. This World class should be intitialized with the maximal growth rate (g), the initial population (pop, and the carrying capacity (pbar). Make sure that it accumulates a history at each step.

class World(object):
    def __init__(self, prms):
        self.g = prms['g']
        self.pop0 = prms['pop']  #initial population
        self.pbar = prms['pbar'] #reference population
    def setup(self):
        self.relpop = self.pop0 / self.pbar
        self.nextLogistic = lambda x: x + self.g * x * (1 - x)
        self.trajectory = [self.pop0]
        self.ticks = 0
    def step(self):
        self.relpop = self.nextLogistic(self.relpop)
        self.ticks += 1
    def recordStep(self):
        newpop = self.pbar * self.relpop  #population
        self.trajectory.append(newpop)    #data collection
    def runsim(self):
        self.setup()
        for _ in range(50):
            self.step()                              #simulation step
            self.trajectory.append(self.pbar * self.relpop) #data collection

Set up and run a simulation as follows.

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

Many Python packages can produce visualizations of the results. For two-dimensional plots, matplotlib is an excellent and popular choice.

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

References

Copyright © 2016–2023 Alan G. Isaac. All rights reserved.

version:

2023-08-09

Copyright © 2016–2023 Alan G. Isaac. All rights reserved.

version:

2023-08-09