# Modeling and Simulation in Python

Chapter 4: Predict

Copyright 2017 Allen Downey

License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)


### Analysis with SymPy


```python
from sympy import *
```

The following line sets up Jupyter notebook to display math.


```python
init_printing() 
```

And this function provides the option of showing results in LaTeX format.


```python
from sympy.printing import latex

def show(expr, show_latex=False):
    if show_latex:
        print(latex(expr))
    return expr
```

Create a symbol for time.


```python
t = symbols('t')
show(t)
```

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


```python
expr = t + 1
show(expr)
```

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

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


```python
f = Function('f')
show(f)
```




    f



The type of `f` is `UndefinedFunction`


```python
type(f)
```




    sympy.core.function.UndefinedFunction



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


```python
f(t)
```

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


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


```python
type(dfdt)
```




    sympy.core.function.Derivative



We need a symbol for `alpha`


```python
alpha = symbols('alpha')
show(alpha)
```

Now we can write the differential equation for proportional growth.


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

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


```python
solution_eq = dsolve(eq1)
show(solution_eq)
```

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


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


```python
particular = solution_eq.subs(C1, p0)
show(particular)
```

In the next example, we'll 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 = symbols('r K')
```

Now we can write the differential equation.


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

And solve it.


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

`rhs` selects the right-hand side of the solution.


```python
general = solution_eq.rhs
show(general)
```

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


```python
at0 = general.subs(t, 0)
show(at0)
```

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

So we'll create the equation `at0 = p0` 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 sequence 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 sequence, so we have to use the bracket operator, `[0]`, to select the first one.


```python
solutions = solve(Eq(at0, p0), C1)
type(solutions), len(solutions)
```




    (list, 1)




```python
value_of_C1 = solutions[0]
show(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 method that tries to simplify it.


```python
particular = simplify(particular)
show(particular)
```

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


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

[In some places](https://en.wikipedia.org/wiki/Population_growth#Logistic_equation) you'll see this solution, which is called the "logistic function" 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 - p0) / p0
show(A)
```


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

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


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

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 is a surprisingly hard problem; in fact, there is no algorithm that can solve it in general.

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


```python
# Solution

solution_eq = dsolve(eq3)
solution_eq
```


```python
# Solution

general = solution_eq.rhs
general
```


```python
# Solution

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


```python
# Solution

solutions = solve(Eq(at0, p0), C1)
value_of_C1 = solutions[0]
value_of_C1
```


```python
# Solution

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

**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 = alpha f(t) (1 - f(t)/K)

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


```python

```


```python

```

cc595bdd7a82a5ec5b6e693fe3715e5f3afab3fd
ipynb
Jupyter Notebook
code/chap04sympy.ipynb
annagriffin/ModSimPy
6f451aae0ce88289734b1153426bb6af933b8dd6
MIT
2021-04-13T01:09:39.000Z
2021-04-13T01:09:39.000Z
code/chap04sympy.ipynb
annagriffin/ModSimPy
6f451aae0ce88289734b1153426bb6af933b8dd6
MIT
code/chap04sympy.ipynb
annagriffin/ModSimPy
6f451aae0ce88289734b1153426bb6af933b8dd6
MIT
2021-04-13T01:10:41.000Z
2021-04-13T01:10:41.000Z
Qwen/Qwen-72B
1. YES
2. YES