Run this code before running the remaining code of the lesson (make sure you have installed "Plots", "LsqFit" and "ODE" packages!!):

In [1]:

```
using Plots;
pyplot()
using LsqFit;
using ODE;
;
```

last update October 10th, 2016

tested ok on Julia 0.4.7

The United Nations considers three possible scenarios of population growth in this century: a. Exponential Growth, b. Capped Growth and c. Declining Population [1]. The growth in each scenario can be forecasted by different variants of a population model called logistic model [2], which is based on ordinary differential equations (ODEs). The logistic model is defined as:

\begin{align*} \frac{dP}{dt} = pgr*P*(1-\frac{P}{K}) \end{align*}where $P$ is the population in billions of people, $t$ is the time in years, $pgr$ is the population growth rate, and $K$ is the carrying capacity in billions of people. In this lesson notebook we show how to use the logistic model to forecast the future world's population under the first two scenarios (a. and b.) using Julia language. In the practice notebook the student will work with the third scenario.

References:

[1] https://en.wikipedia.org/wiki/Projections_of_population_growth

[2] https://en.wikipedia.org/wiki/Logistic_function#In_ecology:_modeling_population_growth

We begin by creating vectors with the world's population data for the periods 1500-2017 and 1900-2017:

In [2]:

```
# historic population data from 1500 to 2017
historic_population = [0.46, 0.58, 0.68, 0.79, 0.98, 1.26, 1.65, 2.52, 5.98, 6.70, 6.90, 7.05, 7.48]
historic_year = [1500, 1600, 1700, 1750, 1800, 1850, 1900, 1950, 1999, 2008, 2010, 2012, 2017]
historic_reduced_population = [1.65, 2.52, 5.98, 6.70, 6.90, 7.05, 7.48]
historic_reduced_year = [1900, 1950, 1999, 2008, 2010, 2012, 2017]
;
```

Reference of world population for 1500 to 1950:

http://www.census.gov/population/international/data/worldpop/table_history.php

Reference of world population for 1950 to 2017:

Now we visualize the population throughout the 1500-2017 period:

In [3]:

```
# Plotting the data
plot(historic_year, historic_population; label="Historic data", marker=([:hex :d]), w=3, background_color=RGB(0.96,0.96,0.96))
xaxis!("Year",(1500,2020),1500:100:2000)
yaxis!("Population (billions)",(0,9),0:2:9)
title!("Historic world population in the period 1500-2017")
```

Out[3]:

We are now ready to begin forecasting the world's population in the coming years!

The Exponential Growth scenario assumes there are so much resources in Earth that the carrying capacity is infinite, $K = \infty$. By substituting $K$ in the logistic model equation we obtain:

\begin{align*} \frac{dP}{dt} = pgr*P \end{align*}The population growth rate $pgr$ is unknown, so we estimate it via least squares using the population data from 1900 to 2017.

In [10]:

```
# Least squares estimation of pgr
historic_reduced_year_shifted = [0, 50, 99, 108, 110, 112, 117]
model(x, p) = 1.65*exp(x.*p[1]); # 1.65 is the population in 1900
fit = curve_fit(model, historic_reduced_year_shifted, historic_reduced_population, [0.01])
pgr = fit.param[1] # estimated population growth rate
;
```

Just to check that the $pgr$ obtained fits well the data, now we integrate the ODE from 1900 to 2017 (using the population of 1900 as the initial population) and then compare results with the historical population during this period.

In [12]:

```
# years used to estimate pgr
year_start = 1900
year_end = 2017
year_diff = year_end - year_start;
```

In [13]:

```
# defining and integrating the ODE
population_change(t,P) = pgr*P[1] # Definition of the ODE
P0 = 1.65 # initial population in billions (corresponding to the world population in 1900)
t = linspace(0,year_diff,year_diff+1) # time span
time, P = ode45(population_change, P0, t;points=:specified); # integrating the ODE
;
```

In [14]:

```
# Plotting the data
plot(historic_year, historic_population; label="Historic data", marker=([:hex :d]), w=3, background_color=RGB(0.96,0.96,0.96))
xaxis!("Year",(1500,2020),1500:100:2000)
yaxis!("Population (billions)",(0,9),0:2:9)
title!("Historic world population in the years 1500-2017")
year_x_axis = linspace(year_start,year_end,year_diff+1)
plot!(year_x_axis, P, label="ODE"; label="ODE", w=3)
```

Out[14]:

Clearly, the $pgr$ obtained fits well the data!

Now we predict the world's population from 2017 to 2100 using the value of $pgr$ found and taking as initial condition the world's population in 2017.

In [15]:

```
# prediction range
prediction_year_start = 2017;
prediction_year_end = 2100;
prediction_year_diff = prediction_year_end - prediction_year_start
;
```

In [16]:

```
# integrating the ODE
P_0 = 7.48; # initial population in billions (corresponding to the world's population in 2017)
t = linspace(0, prediction_year_diff, prediction_year_diff+1) # time span
time, P = ode45(population_change, P_0, t;points=:specified);
```

In [17]:

```
# Visualizing historic and predicted data
plot(historic_reduced_year, historic_reduced_population; label="Historic data", marker=([:hex :d]), w=3, background_color=RGB(0.96,0.96,0.96))
xaxis!("Year",(1900,2100),1900:50:2100)
yaxis!("Population (billions)",(0,25),0:5:25)
title!("Predicted world population")
prediction_year_x = linspace(prediction_year_start, prediction_year_end, prediction_year_diff+1)
plot!(prediction_year_x, P, label="ODE"; label="ODE", w=3)
```

Out[17]:

The model predicts that there will be near 22 billion people in the year 2100! The predicted population is very high because the model does not consider resource scarcity.

Since resources available on Earth are finite, the world's population in reality would not sustain exponential growth indefinitely. The following variant of the logistic model takes into account this, and assumes the carrying capacity (denoted by $K$) is an upper population threshold (i.e. population cannot grow bigger than it). The ODE is:

\begin{align*} \frac{dP}{dt} = pgr*P*(1-\frac{P}{K}) \end{align*}where $P$ is the population in billions of people, $t$ is the time in years, $pgr$ is the population growth rate, and $K$ is the carrying capacity in billions of people. In the following, we predict the world's population in the period 2017-2100 using this model. We use the value of $pgr$ found in last section and take $K$ to be equal to 11 billion people.

In [18]:

```
# Parameters of the ODE
pgr = 0.013
K = 11 # population cap in billions of people (carrying capacity)
# ODE definition
population_change_with_cap(t,P) = pgr*P*(1-P/K)
# Initial condition
P_0 = 7.48 # world's population in 2017
# Time span
t = linspace(0, prediction_year_diff, prediction_year_diff+1)
# integrating the ODE
time, P = ode45(population_change_with_cap, P_0, t;points=:specified);
;
```

In [19]:

```
# Visualizing historic and predicted data
# Plotting the data
plot(historic_reduced_year, historic_reduced_population; label="Historic data", marker=([:hex :d]), w=3, background_color=RGB(0.96,0.96,0.96))
xaxis!("Year",(1900,2100),1900:50:2100)
yaxis!("Population (billions)",(0,12),0:2:12)
title!("Historic and predicted world population")
prediction_year_x = linspace(prediction_year_start, prediction_year_end, prediction_year_diff+1)
plot!(prediction_year_x, P, label="ODE", w=3)
```

Out[19]:

The predicted population for the year 2100 is near 9.5 billion people, and in future centuries it will end up converging to 11 billion (carrying capacity).