# Chapter 9

*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
# check if the libraries we need are installed

try:
    import pint
except ImportError:
    !pip install pint
    
try:
    from modsim import *
except ImportError:
    !pip install modsimpy
```

In this chapter we express the models from previous chapters as
difference equations and differential equations, solve the equations,
and derive the functional forms of the solutions. We also discuss 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_0$ is a given initial population, and $c$ is constant annual growth.
This way of representing the model 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 it is not hard to conclude 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$$ 

This result is 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 model is **discrete**, which
means it 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 a **continuous** model might be more realistic.

We can make this model continuous by writing the rate of change in the
form of a derivative: 

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

This way of representing the model is a **differential equation**; see <http://modsimpy.com/diffeq>.

We can solve this differential equation if 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; 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 growth rate versus population and saw net crossed 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 the
    analysis in Section xxx, we get the general result that
    $K=-\alpha/\beta$.

-   Analysis often 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 type

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

Before we can use SymPy, we have to import it:

SymPy defines a `Symbol` object that represents symbolic variable names,
functions, and other mathematical entities.

The `symbols` function takes a string and returns `Symbol` objects. So
if we run this assignment:


```python
from sympy import symbols

t = symbols('t')
```

Python understands that `t` is a symbol, not a numerical value. If we
now run


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




$\displaystyle t + 1$



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




$\displaystyle 3$



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


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


```python
f(t)
```




$\displaystyle f{\left(t \right)}$



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




$\displaystyle \frac{d}{d t} f{\left(t \right)}$



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




$\displaystyle \frac{d}{d t} f{\left(t \right)} = \alpha f{\left(t \right)}$



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




$\displaystyle f{\left(t \right)} = C_{1} e^{\alpha t}$



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




$\displaystyle f{\left(t \right)} = p_{0} e^{\alpha t}$



The result is called the **exponential growth curve**; see
<http://modsimpy.com/expo>.

## Solving the quadratic growth model

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




$\displaystyle \frac{d}{d t} f{\left(t \right)} = r \left(1 - \frac{f{\left(t \right)}}{K}\right) f{\left(t \right)}$



And solve it.


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




$\displaystyle f{\left(t \right)} = \frac{K e^{C_{1} K + r t}}{e^{C_{1} K + r t} - 1}$



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


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




$\displaystyle \frac{K e^{C_{1} K + r t}}{e^{C_{1} K + r t} - 1}$



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


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




$\displaystyle \frac{K e^{C_{1} K}}{e^{C_{1} K} - 1}$



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.  In this case, [we have reason to expect only one solution](https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem), but we still get a list, so we have to use the bracket operator, `[0]`, to select the first one.


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




    (list, 1)




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




$\displaystyle \frac{\log{\left(- \frac{p_{0}}{K - p_{0}} \right)}}{K}$



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




$\displaystyle - \frac{K p_{0} e^{r t}}{\left(K - p_{0}\right) \left(- \frac{p_{0} e^{r t}}{K - p_{0}} - 1\right)}$



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


```python
particular.simplify()
```




$\displaystyle \frac{K p_{0} e^{r t}}{K + p_{0} e^{r t} - p_{0}}$



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




$\displaystyle \frac{K - p_{0}}{p_{0}}$




```python
from sympy import exp

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




$\displaystyle \frac{K}{1 + \frac{\left(K - p_{0}\right) e^{- r t}}{p_{0}}}$



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


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




$\displaystyle 0$



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

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 + \alpha t$$

where $\alpha$ is the growth rate.

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

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

Finally, the quadratic growth model is called **logistic growth** because the solution is the 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.

## Exercises

**Exercise:** Solve the quadratic growth equation using the alternative parameterization

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


```python
# Solution

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


```python
# Solution

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




$\displaystyle \frac{d}{d t} f{\left(t \right)} = \alpha f{\left(t \right)} + \beta f^{2}{\left(t \right)}$




```python
# Solution

solution_eq = dsolve(eq3)
solution_eq
```




$\displaystyle f{\left(t \right)} = \frac{\alpha e^{\alpha \left(C_{1} + t\right)}}{\beta \left(1 - e^{\alpha \left(C_{1} + t\right)}\right)}$




```python
# Solution

general = solution_eq.rhs
general
```




$\displaystyle \frac{\alpha e^{\alpha \left(C_{1} + t\right)}}{\beta \left(1 - e^{\alpha \left(C_{1} + t\right)}\right)}$




```python
# Solution

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


```python
# Solution

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




$\displaystyle \frac{\log{\left(\frac{\beta p_{0}}{\alpha + \beta p_{0}} \right)}}{\alpha}$




```python
# Solution

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




$\displaystyle \frac{\alpha p_{0} e^{\alpha t}}{\alpha - \beta p_{0} e^{\alpha t} + \beta p_{0}}$



**Exercise:**  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`.

87d7b420dfb6e66ab71cf6a0a36b854d9eae2df7
ipynb
Jupyter Notebook
jupyter/chap09.ipynb
SadPuffin/ModSimPy
6cac03b114a0bf0893aa160048ba6dae460f5f86
MIT
jupyter/chap09.ipynb
SadPuffin/ModSimPy
6cac03b114a0bf0893aa160048ba6dae460f5f86
MIT
jupyter/chap09.ipynb
SadPuffin/ModSimPy
6cac03b114a0bf0893aa160048ba6dae460f5f86
MIT
Qwen/Qwen-72B
1. YES
2. YES