# Analysis of Population Growth

*Modeling and Simulation in Python*

Copyright 2021 Allen Downey

License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)


```python
# install Pint if necessary

try:
    import pint
except ImportError:
    !pip install pint
```


```python
# download modsim.py if necessary

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
    
download('https://raw.githubusercontent.com/AllenDowney/' +
         'ModSimPy/master/modsim.py')
```


```python
# import functions from modsim

from modsim import *
```

In this chapter we'll express the models from previous chapters as
difference equations and differential equations, solve the equations,
and derive the functional forms of the solutions. And I'll present some thoughts about the complementary roles of mathematical analysis and simulation.

## Recurrence Relations

The population models in the previous chapter and this one are simple
enough that we didn't really need to run simulations. We could have
solved them mathematically. For example, we wrote the constant growth
model like this:

```
results[t+1] = results[t] + annual_growth
```

In mathematical notation, we would write the same model like this:

$$x_{n+1} = x_n + c$$ 

where $x_n$ is the population during year $n$, 
$x_{n+1}$ is the population during year $n+1$,
and $c$ is constant annual growth.
This way of representing the model, where future population is a function of current population, is a *recurrence relation*; see
<http://modsimpy.com/recur>.

Sometimes it is possible to solve a recurrence relation by writing an
equation that computes $x_n$, for a given value of $n$, directly; that
is, without computing the intervening values from $x_1$ through
$x_{n-1}$.

In the case of constant growth we can see that $x_1 = x_0 + c$, and
$x_2 = x_1 + c$. Combining these, we get $x_2 = x_0 + 2c$, then
$x_3 = x_0 + 3c$, and we can see that in general

$$x_n = x_0 + nc$$ 

So if we want to know $x_{100}$ and we don't care
about the other values, we can compute it with one multiplication and
one addition.

We can also write the proportional model as a recurrence relation:

$$x_{n+1} = x_n + \alpha x_n$$ 

Or more conventionally as:

$$x_{n+1} = x_n (1 + \alpha)$$ 

Now we can see that $x_1 = x_0 (1 + \alpha)$, and $x_2 = x_0 (1 + \alpha)^2$, and in general

$$x_n = x_0 (1 + \alpha)^n$$

A sequence with this functional form is called a *geometric progression*;
see <http://modsimpy.com/geom>. When $\alpha$ is positive, the factor
$1+\alpha$ is greater than 1, so the elements of the sequence grow
without bound.

Finally, we can write the quadratic model like this:

$$x_{n+1} = x_n + \alpha x_n + \beta x_n^2$$ 

or with the more conventional parameterization like this:

$$x_{n+1} = x_n + r x_n (1 - x_n / K)$$ 

There is no analytic solution to this equation, but we can approximate it with a differential equation and solve that, which is what we'll do in the next section.

## Differential Equations

Starting again with the constant growth model 

$$x_{n+1} = x_n + c$$ 

If we define $\Delta x$ to be the change in $x$ from one time step to the next, we can write: 

$$\Delta x = x_{n+1} - x_n = c$$ 

If we define
$\Delta t$ to be the time step, which is one year in the example, we can
write the rate of change per unit of time like this:

$$\frac{\Delta x}{\Delta t} = c$$ 

This is a *discrete* model, which
means time is only defined at integer values of $n$ and not in between.
But in reality, people are born and die all the time, not once a year,
so it might be more realistic use a *continuous* model, which means
time is defined at all values of $t$, not just integers. 

In a continuous model, we write the rate of change in the
form of a derivative: 

$$\frac{dx}{dt} = c$$ 

This way of representing the model is a *differential equation*, which is an equation that involves at least one derivative (see <http://modsimpy.com/diffeq>).
To solve this equation, we multiply both sides by $dt$: 

$$dx = c dt$$ 

And then integrate both sides: 

$$x(t) = c t + x_0$$

Similarly, we can write the proportional growth model like this:

$$\frac{\Delta x}{\Delta t} = \alpha x$$ 

And as a differential equation
like this: 

$$\frac{dx}{dt} = \alpha x$$ 

If we multiply both sides by
$dt$ and divide by $x$, we get 

$$\frac{1}{x}~dx = \alpha~dt$$ 

Now we
integrate both sides, yielding: 

$$\ln x = \alpha t + K$$ 

where $\ln$ is the natural logarithm and $K$ is the constant of integration.

Exponentiating both sides, we have

$$\exp(\ln(x)) = \exp(\alpha t + K)$$ 

The exponential function can be written $\exp(x)$ or $e^x$. In this book I use the first form because it resembles the Python code.
We can rewrite the previous equation as

$$x = \exp(\alpha t) \exp(K)$$ 

Since $K$ is an arbitrary constant,
$\exp(K)$ is also an arbitrary constant, so we can write

$$x = C \exp(\alpha t)$$ 

where $C = \exp(K)$. There are many solutions to this differential equation, with different values of $C$. The particular solution we want is the one that has the value $x_0$ when $t=0$.

When $t=0$, $x(t) = C$, so $C = x_0$ and the solution we want is

$$x(t) = x_0 \exp(\alpha t)$$ 

If you would like to see this derivation done more carefully, you might like this video:
<http://modsimpy.com/khan1>.

## Analysis and Simulation

Once you have designed a model, there are generally two ways to proceed: simulation and analysis. Simulation often comes in the form of a computer program that models changes in a system over time, like births and deaths, or bikes moving from place to place. Analysis often comes in the form of algebra and calculus; that is, symbolic manipulation using mathematical notation.

Analysis and simulation have different capabilities and limitations.
Simulation is generally more versatile; it is easy to add and remove parts of a program and test many versions of a model, as we have done in the previous examples.

But there are several things we can do with analysis that are harder or impossible with simulations:

-   With analysis we can sometimes compute, exactly and efficiently, a
    value that we could only approximate, less efficiently, with
    simulation. For example, in the quadratic model we plotted net growth versus population and saw it crosses through zero when the population is
    near 14 billion.  We could estimate the crossing point using a
    numerical search algorithm (more about that later). But with a bit of algebra, we derived the general result that $K=-\alpha/\beta$.

-   Analysis sometimes provides "computational shortcuts", that is, the
    ability to jump forward in time to compute the state of a system
    many time steps in the future without computing the intervening
    states.

-   We can use analysis to state and prove generalizations about models;
    for example, we might prove that certain results will always or
    never occur. With simulations, we can show examples and sometimes
    find counterexamples, but it is hard to write proofs.

-   Analysis can provide insight into models and the systems they
    describe; for example, sometimes we can identify regimes of
    qualitatively different behavior and key parameters that control
    those behaviors.

When people see what analysis can do, they sometimes get drunk with
power and imagine that it gives them a special ability to see past the veil of the material world and discern the laws of mathematics that govern the universe. When they analyze a model of a physical system, they talk about "the math behind it" as if our world is the mere shadow of a world of ideal mathematical entities (I am not making this up; see <http://modsimpy.com/plato>.).

This is, of course, nonsense. Mathematical notation is a language
designed by humans for a purpose, specifically to facilitate symbolic
manipulations like algebra. Similarly, programming languages are
designed for a purpose, specifically to represent computational ideas
and run programs.

Each of these languages is good for the purposes it was designed for and less good for other purposes. But they are often complementary, and one of the goals of this book is to show how they can be used together.

## Analysis with WolframAlpha

Until recently, most analysis was done by rubbing graphite on wood
pulp, a process that is laborious and error-prone. A useful
alternative is symbolic computation. If you have used services like
WolframAlpha, you have used symbolic computation.

For example, if you go to <https://www.wolframalpha.com/> and enter

```
df(t) / dt = alpha f(t)
```

WolframAlpha infers that `f(t)` is a function of `t` and `alpha` is a
parameter; it classifies the query as a "first-order linear ordinary
differential equation", and reports the general solution:

$$f(t) = c_1 \exp(\alpha t)$$ 

If you add a second equation to specify the initial condition:

```
df(t) / dt = alpha f(t),  f(0) = p_0
```

WolframAlpha reports the particular solution:

$$f(t) = p_0 \exp(\alpha t)$$

WolframAlpha is based on Mathematica, a powerful programming language
designed specifically for symbolic computation.

## Analysis with SymPy

Python has a library called SymPy that provides symbolic computation
tools similar to Mathematica. They are not as easy to use as
WolframAlpha, but they have some other advantages.

To use it, we'll define `Symbol` objects that represent variable names and functions.
The `symbols` function takes a string and returns `Symbol` objects.
So if we run this assignment:


```python
from sympy import symbols

t = symbols('t')
```

Now when we use `t`, Python understands that it is a symbol, not a numerical value. For example, if we use `t` as part of an expression, like this,


```python
expr = t + 1
expr
```

Python doesn't try to perform numerical addition; rather, it creates a
new `Symbol` that represents the sum of `t` and `1`. We can evaluate the
sum using `subs`, which substitutes a value for a symbol. This example
substitutes 2 for `t`:


```python
expr.subs(t, 2)
```

Functions in SymPy are represented by a special kind of `Symbol`:


```python
from sympy import Function

f = Function('f')
f
```

Now if we write `f(t)`, we get an object that represents the evaluation of a function, $f$, at a value, $t$. 


```python
f(t)
```

But again SymPy doesn't actually
try to evaluate it.

## Differential Equations In SymPy

SymPy provides a function, `diff`, that can differentiate a function. We can apply it to `f(t)` like this:


```python
from sympy import diff

dfdt = diff(f(t), t)
dfdt
```

The result is a `Symbol` that represents the derivative of `f` with
respect to `t`. But again, SymPy doesn't try to compute the derivative
yet.

To represent a differential equation, we use `Eq`:


```python
from sympy import Eq

alpha = symbols('alpha')
eq1 = Eq(dfdt, alpha*f(t))
eq1
```

The result is an object that represents an equation.  Now
we can use `dsolve` to solve this differential equation:


```python
from sympy import dsolve

solution_eq = dsolve(eq1)
solution_eq
```

The result is the *general
solution*, which still contains an unspecified constant, $C_1$. To get the *particular solution* where $f(0) = p_0$, we substitute `p_0` for `C1`. First, we have to create two more symbols:


```python
C1, p_0 = symbols('C1 p_0')
```

Now we can perform the substitution:


```python
particular = solution_eq.subs(C1, p_0)
particular
```

When $t=0$, $e^{\alpha t}$ is $1$, so $f(t)$ is $p_0$, which is what we wanted.

## Solving the Quadratic Growth Model

To solve the quadratic growth curve, we'll use the `r, K` parameterization, so we'll need two more symbols:


```python
r, K = symbols('r K')
```

Now we can write the differential equation.


```python
eq2 = Eq(diff(f(t), t), r * f(t) * (1 - f(t)/K))
eq2
```

And solve it.


```python
solution_eq = dsolve(eq2)
solution_eq
```

The result, `solution_eq`, contains `rhs`, which is the right-hand side of the solution.


```python
general = solution_eq.rhs
general
```

We can evaluate the right-hand side at $t=0$


```python
at_0 = general.subs(t, 0)
at_0
```

Now we want to find the value of `C1` that makes `f(0) = p_0`.

So we'll create the equation `at_0 = p_0` and solve for `C1`.  Because this is just an algebraic identity, not a differential equation, we use `solve`, not `dsolve`.


```python
from sympy import solve

solutions = solve(Eq(at_0, p_0), C1)
```

The result from `solve` is a list of solutions.  


```python
type(solutions), len(solutions)
```

In this case, there is only one solution, but we still get a list, so we have to use the bracket operator, `[0]`, to select the first one.


```python
value_of_C1 = solutions[0]
value_of_C1
```

Now in the general solution, we want to replace `C1` with the value of `C1` we just figured out.


```python
particular = general.subs(C1, value_of_C1)
particular
```

The result is complicated, but SymPy provides a function that tries to simplify it.


```python
simpler = particular.simplify()
simpler
```

This function is called the *logistic growth curve*; see
<http://modsimpy.com/logistic>. In the context of growth models, the
logistic function is often written like this:

$$f(t) = \frac{K}{1 + A \exp(-rt)}$$ 

where $A = (K - p_0) / p_0$.

We can use SymPy to confirm that these two forms are equivalent.  First we represent the alternative version of the logistic function:


```python
A = (K - p_0) / p_0
A
```


```python
from sympy import exp

logistic = K / (1 + A * exp(-r*t))
logistic
```

To see whether two expressions are equivalent, we can check whether their difference simplifies to 0.


```python
(particular - logistic).simplify()
```

This test only works one way: if SymPy says the difference reduces to 0, the expressions are definitely equivalent (and not just numerically close).

But if SymPy can't find a way to simplify the result to 0, that doesn't necessarily mean there isn't one.  Testing whether two expressions are equivalent is a surprisingly hard problem; in fact, there is no algorithm that can solve it in general.

If you use SymPy to compute and expression, and then want to evaluate that expression in Python, SumPy provides a function called `pycode` that generates Python code:


```python
from sympy.printing.pycode import pycode

pycode(simpler)
```

If you would like to see this differential equation solved by hand, you might like this video: <http://modsimpy.com/khan2>

## Summary

In this chapter we wrote the growth models from the previous chapters in terms of difference and differential equations. We solved some of these equations by hand; for others, we used WolframAlpha and SymPy.

What I called the "constant growth" model is more commonly called *linear growth* because the solution is a line.  If we model time as continuous, the solution is

$$f(t) = p_0 + c t$$

where $c$ is net annual growth.

Similarly, the proportional growth model is usually called *exponential growth* because the solution is an exponential function:

$$f(t) = p_0 \exp{\alpha t}$$

Finally, the quadratic growth model is called *logistic growth* because the solution is a logistic function:

$$f(t) = \frac{K}{1 + A \exp(-rt)}$$ 

where $A = (K - p_0) / p_0$.

I avoided these terms until now because they are based on results we had not derived yet.

With that, we are done modeling world population growth.
The next chapter presents case studies where you can apply the tools we have learned so far.

## Exercises

This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. 
You can access the notebooks at <https://allendowney.github.io/ModSimPy/>.

### Exercise 1

 Use SymPy to solve the quadratic growth equation using the alternative parameterization

$$ \frac{df(t)}{dt} = \alpha f(t) + \beta f^2(t) $$



```python
# Solution goes here
```


```python
# Solution goes here
```


```python
# Solution goes here
```


```python
# Solution goes here
```


```python
# Solution goes here
```


```python
# Solution goes here
```


```python
# Solution goes here
```

### Exercise 2

  Use [WolframAlpha](https://www.wolframalpha.com/) to solve the quadratic growth model, using either or both forms of parameterization:

    df(t) / dt = alpha f(t) + beta f(t)^2

or

    df(t) / dt = r f(t) (1 - f(t)/K)

Find the general solution and also the particular solution where `f(0) = p_0`.


```python

```

6a4b6a95d04b10946d5d7006c2db685e111d0b0f
ipynb
Jupyter Notebook
chapters/chap09.ipynb
ABajad/ModSimPy
cc14dbb22c8c11ab3f25477d04575709a955812a
MIT
chapters/chap09.ipynb
ABajad/ModSimPy
cc14dbb22c8c11ab3f25477d04575709a955812a
MIT
chapters/chap09.ipynb
ABajad/ModSimPy
cc14dbb22c8c11ab3f25477d04575709a955812a
MIT
Qwen/Qwen-72B
1. YES
2. YES