<jupyter_start><jupyter_text>Chapter 9 *Modeling and Simulation in Python*Copyright 2021 <NAME>License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)<jupyter_code># 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 *<jupyter_output><empty_output><jupyter_text>The following displays SymPy expressions and provides the option of showing results in LaTeX format.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Analysis with SymPy Create a symbol for time.<jupyter_code>import sympy as sp

t = sp.symbols('t')
t<jupyter_output><empty_output><jupyter_text>If you combine symbols and numbers, you get symbolic expressions.<jupyter_code>expr = t + 1
expr<jupyter_output><empty_output><jupyter_text>The result is an `Add` object, which just represents the sum without trying to compute it.<jupyter_code>type(expr)<jupyter_output><empty_output><jupyter_text>`subs` can be used to replace a symbol with a number, which allows the addition to proceed.<jupyter_code>expr.subs(t, 2)<jupyter_output><empty_output><jupyter_text>`f` is a special class of symbol that represents a function.<jupyter_code>f = sp.Function('f')
f<jupyter_output><empty_output><jupyter_text>The type of `f` is `UndefinedFunction`<jupyter_code>type(f)<jupyter_output><empty_output><jupyter_text>SymPy understands that `f(t)` means `f` evaluated at `t`, but it doesn't try to evaluate it yet.<jupyter_code>f(t)<jupyter_output><empty_output><jupyter_text>`diff` returns a `Derivative` object that represents the time derivative of `f`<jupyter_code>dfdt = sp.diff(f(t), t)
dfdt
type(dfdt)<jupyter_output><empty_output><jupyter_text>We need a symbol for `alpha`<jupyter_code>alpha = sp.symbols('alpha')
alpha<jupyter_output><empty_output><jupyter_text>Now we can write the differential equation for proportional growth.<jupyter_code>eq1 = sp.Eq(dfdt, alpha*f(t))
eq1<jupyter_output><empty_output><jupyter_text>And use `dsolve` to solve it.  The result is the general solution.<jupyter_code>solution_eq = sp.dsolve(eq1)
solution_eq<jupyter_output><empty_output><jupyter_text>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`<jupyter_code>C1, p_0 = sp.symbols('C1 p_0')
particular = solution_eq.subs(C1, p_0)
particular<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>r, K = sp.symbols('r K')<jupyter_output><empty_output><jupyter_text>Now we can write the differential equation.<jupyter_code>eq2 = sp.Eq(sp.diff(f(t), t), r * f(t) * (1 - f(t)/K))
eq2<jupyter_output><empty_output><jupyter_text>And solve it.<jupyter_code>solution_eq = sp.dsolve(eq2)
solution_eq<jupyter_output><empty_output><jupyter_text>The result, `solution_eq`, contains `rhs`, which is the right-hand side of the solution.<jupyter_code>general = solution_eq.rhs
general<jupyter_output><empty_output><jupyter_text>We can evaluate the right-hand side at $t=0$<jupyter_code>at_0 = general.subs(t, 0)
at_0<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>solutions = sp.solve(sp.Eq(at_0, p_0), C1)
type(solutions), len(solutions)
value_of_C1 = solutions[0]
value_of_C1<jupyter_output><empty_output><jupyter_text>Now in the general solution, we want to replace `C1` with the value of `C1` we just figured out.<jupyter_code>particular = general.subs(C1, value_of_C1)
particular<jupyter_output><empty_output><jupyter_text>The result is complicated, but SymPy provides a method that tries to simplify it.<jupyter_code>particular = sp.simplify(particular)
particular<jupyter_output><empty_output><jupyter_text>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`<jupyter_code>particular.subs(t, 0)<jupyter_output><empty_output><jupyter_text>This solution is called the [logistic function](https://en.wikipedia.org/wiki/Population_growthLogistic_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:<jupyter_code>A = (K - p_0) / p_0
A
logistic = K / (1 + A * sp.exp(-r*t))
logistic<jupyter_output><empty_output><jupyter_text>To see whether two expressions are equivalent, we can check whether their difference simplifies to 0.<jupyter_code>sp.simplify(particular - logistic)<jupyter_output><empty_output><jupyter_text>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) $<jupyter_code>beta = sp.symbols('beta')
beta
eq3 = sp.Eq(sp.diff(f(t), t), alpha * f(t) + beta * f(t)**2)
eq3
solution_eq3 = sp.dsolve(eq3)
solution_eq3
general3 = solution_eq3.rhs
general3
at_03 = general3.subs(t, 0)
at_03
solutions3 = sp.solve(sp.Eq(at_03, p_0), C1)
solutions3[0]
particular3 = sp.simplify(general3.subs(C1, solutions3[0]))
particular3<jupyter_output><empty_output><jupyter_text>**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)^2or    df(t) / dt = r f(t) (1 - f(t)/K)Find the general solution and also the particular solution where `f(0) = p_0`.<jupyter_code>Please see the solution above<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10679196924391286, 'fraction_numerical': 0.014238929232521715, 'mean_word_length': 6.339602925809823, 'pattern_counts': {'":': 0, '<': 110, '<?xml version=': 0, '>': 110, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '9795', 'license': 'MIT', 'n_tokens_mistral': 2472, 'n_tokens_neox': 2385, 'n_words': 700, 'path': 'colab/chap09.ipynb', 'repo_name': 'ferrysany/ModSimPy', 'size': '7029'}