# Analysis of Population Growth#

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.

## Difference Equations#

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:

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 *difference equation*; see
https://en.wikipedia.org/wiki/Linear_difference_equation.

For a given value of \(n\), sometimes it is possible to compute \(x_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

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 difference equation:

Or more conventionally as:

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

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:

or with the more conventional parameterization like this:

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

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

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:

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 to 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:

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

And then integrate both sides:

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

And as a differential equation like this:

If we multiply both sides by \(dt\) and divide by \(x\), we get

Now we integrate both sides, yielding:

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

Exponentiating both sides, we have

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

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

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

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 qualitatively different ways the system can behave 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:

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:

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 names of variables and functions.
The `symbols`

function takes a string and returns `Symbol`

objects.

```
from sympy import symbols
t = symbols('t')
```

Now when we use `t`

, Python treats it like a variable name rather than a specific number.
For example, if we use `t`

as part of an expression, like this,

```
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`

:

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

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

:

```
from sympy import Function
f = Function('f')
f
```

```
f
```

Now if we write `f(t)`

, we get an object that represents the evaluation of a function, \(f\), at a value, \(t\).

```
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:

```
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`

:

```
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:

```
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 tell Python that `C1`

is a symbol.

```
C1 = symbols('C1')
```

Now we can substitute the value of \(p_0\) for `C1`

.
For example, if \(p_0\) is 1000:

```
particular = solution_eq.subs(C1, 1000)
particular
```

When \(t=0\), the value of \(f(0)\) is \(p_0\), which confirms that this is the solution we want.

```
particular.subs(t, 0)
```

## 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:

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

Now we can write the differential equation.

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

And solve it.

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

The result, `solution_eq`

, contains `rhs`

, which is the right-hand side of the solution.

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

We can evaluate the right-hand side at \(t=0\)

```
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 equation, not a differential equation, we use `solve`

, not `dsolve`

.

```
from sympy import solve
p_0 = symbols('p_0')
solutions = solve(Eq(at_0, p_0), C1)
```

The result from `solve`

is a list of solutions.

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

```
(list, 1)
```

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.

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

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

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

```
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:

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

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

where \(c\) is net annual growth.

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

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

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

## Show code cell content

```
# Solution
alpha, beta = symbols('alpha beta')
```

## Show code cell content

```
# Solution
eq3 = Eq(diff(f(t), t), alpha*f(t) + beta*f(t)**2)
eq3
```

## Show code cell content

```
# Solution
solution_eq = dsolve(eq3)
solution_eq
```

## Show code cell content

```
# Solution
general = solution_eq.rhs
general
```

## Show code cell content

```
# Solution
at_0 = general.subs(t, 0)
```

## Show code cell content

```
# Solution
solutions = solve(Eq(at_0, p_0), C1)
value_of_C1 = solutions[0]
value_of_C1
```

## Show code cell content

```
# Solution
particular = general.subs(C1, value_of_C1)
particular.simplify()
```

### Exercise 2#

Use WolframAlpha 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`

.