# 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
    import pint
    
try:
    from modsim import *
except ImportError:
    !pip install modsimpy
    from modsim import *
```

The following displays SymPy expressions and provides the option of showing results in LaTeX format.


```python
from sympy.printing import latex

def show(expr, show_latex=False):
    """Display a SymPy expression.
    
    expr: SymPy expression
    show_latex: boolean
    """
    if show_latex:
        print(latex(expr))
    return expr
```

### Analysis with SymPy

Create a symbol for time.


```python
import sympy as sp

t = sp.symbols('t')
t
```




$\displaystyle t$



If you combine symbols and numbers, you get symbolic expressions.


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




$\displaystyle t + 1$



The result is an `Add` object, which just represents the sum without trying to compute it.


```python
type(expr)
```




    sympy.core.add.Add



`subs` can be used to replace a symbol with a number, which allows the addition to proceed.


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




$\displaystyle 3$



`f` is a special class of symbol that represents a function.


```python
f = sp.Function('f')
f
```




    f



The type of `f` is `UndefinedFunction`


```python
type(f)
```




    sympy.core.function.UndefinedFunction



SymPy understands that `f(t)` means `f` evaluated at `t`, but it doesn't try to evaluate it yet.


```python
f(t)
```




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



`diff` returns a `Derivative` object that represents the time derivative of `f`


```python
dfdt = sp.diff(f(t), t)
dfdt
```




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




```python
type(dfdt)
```




    sympy.core.function.Derivative



We need a symbol for `alpha`


```python
alpha = sp.symbols('alpha')
alpha
```




$\displaystyle \alpha$



Now we can write the differential equation for proportional growth.


```python
eq1 = sp.Eq(dfdt, alpha*f(t))
eq1
```




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



And use `dsolve` to solve it.  The result is the general solution.


```python
solution_eq = sp.dsolve(eq1)
solution_eq
```




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



We can tell it's a general solution because it contains an unspecified constant, `C1`.

In this example, finding the particular solution is easy: we just replace `C1` with `p_0`


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


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




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



In the next example, we have to work a little harder to find the particular solution.

### Solving the quadratic growth equation 

We'll use the (r, K) parameterization, so we'll need two more symbols:


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

Now we can write the differential equation.


```python
eq2 = sp.Eq(sp.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 = sp.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`.

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
solutions = sp.solve(sp.Eq(at_0, p_0), C1)
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 method that tries to simplify it.


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




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



Often simplicity is in the eye of the beholder, but that's about as simple as this expression gets.

Just to double-check, we can evaluate it at `t=0` and confirm that we get `p_0`


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




$\displaystyle p_{0}$



This solution is called the [logistic function](https://en.wikipedia.org/wiki/Population_growth#Logistic_equation).

In some places you'll see it written in a different form:

$f(t) = \frac{K}{1 + A e^{-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
logistic = K / (1 + A * sp.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
sp.simplify(particular - logistic)
```




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

### Exercises

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

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


```python
beta = sp.symbols('beta')
beta
```




$\displaystyle \beta$




```python
eq3 = sp.Eq(sp.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_eq3 = sp.dsolve(eq3)
solution_eq3
```




$\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
general3 = solution_eq3.rhs
general3
```




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




```python
at_03 = general3.subs(t, 0)
at_03
```




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




```python
solutions3 = sp.solve(sp.Eq(at_03, p_0), C1)
solutions3[0]
```




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




```python
particular3 = sp.simplify(general3.subs(C1, solutions3[0]))
particular3
```




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


```python
Please see the solution above
```

4a83b9f6182251234e0e64421bbd5dd9a42583ab
ipynb
Jupyter Notebook
colab/chap09.ipynb
ferrysany/ModSimPy
4b68634ee847102ad3b0f1816ac0b1c3125018d3
MIT
colab/chap09.ipynb
ferrysany/ModSimPy
4b68634ee847102ad3b0f1816ac0b1c3125018d3
MIT
colab/chap09.ipynb
ferrysany/ModSimPy
4b68634ee847102ad3b0f1816ac0b1c3125018d3
MIT
Qwen/Qwen-72B
1. YES
2. YES