Modeling and Simulation in Python

Chapter 9

License: Creative Commons Attribution 4.0 International



In [1]:

    
# import everything from SymPy.
from sympy import *

# Set up Jupyter notebook to display math.
init_printing()

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



In [2]:

    
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.



In [3]:

    
t = symbols('t')
t









    Out[3]:





$\displaystyle t$

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



In [4]:

    
expr = t + 1
expr









    Out[4]:





$\displaystyle t + 1$

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



In [5]:

    
type(expr)









    Out[5]:





sympy.core.add.Add

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



In [6]:

    
expr.subs(t, 2)









    Out[6]:





$\displaystyle 3$

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



In [7]:

    
f = Function('f')
f









    Out[7]:





f

The type of f is UndefinedFunction



In [8]:

    
type(f)









    Out[8]:





sympy.core.function.UndefinedFunction

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



In [9]:

    
f(t)









    Out[9]:





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

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



In [10]:

    
dfdt = diff(f(t), t)
dfdt









    Out[10]:





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



In [11]:

    
type(dfdt)









    Out[11]:





sympy.core.function.Derivative

We need a symbol for alpha



In [12]:

    
alpha = symbols('alpha')
alpha









    Out[12]:





$\displaystyle \alpha$

Now we can write the differential equation for proportional growth.



In [13]:

    
eq1 = Eq(dfdt, alpha*f(t))
eq1









    Out[13]:





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



In [14]:

    
solution_eq = dsolve(eq1)
solution_eq









    Out[14]:





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



In [15]:

    
C1, p_0 = symbols('C1 p_0')



In [16]:

    
particular = solution_eq.subs(C1, p_0)
particular









    Out[16]:





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



In [17]:

    
r, K = symbols('r K')

Now we can write the differential equation.



In [18]:

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









    Out[18]:





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



In [19]:

    
solution_eq = dsolve(eq2)
solution_eq









    Out[19]:





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



In [20]:

    
general = solution_eq.rhs
general









    Out[20]:





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



In [21]:

    
at_0 = general.subs(t, 0)
at_0









    Out[21]:





$\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, but we still get a list, so we have to use the bracket operator, [0], to select the first one.



In [22]:

    
solutions = solve(Eq(at_0, p_0), C1)
type(solutions), len(solutions)









    Out[22]:





(list, 1)



In [23]:

    
value_of_C1 = solutions[0]
value_of_C1









    Out[23]:





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



In [24]:

    
particular = general.subs(C1, value_of_C1)
particular









    Out[24]:





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



In [25]:

    
particular = simplify(particular)
particular









    Out[25]:





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



In [26]:

    
particular.subs(t, 0)









    Out[26]:





$\displaystyle p_{0}$

This solution is called the logistic function.

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:



In [27]:

    
A = (K - p_0) / p_0
A









    Out[27]:





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



In [28]:

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









    Out[28]:





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



In [29]:

    
simplify(particular - logistic)









    Out[29]:





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



In [30]:

    
# Solution

alpha, beta = symbols('alpha beta')



In [31]:

    
# Solution

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









    Out[31]:





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



In [32]:

    
# Solution

solution_eq = dsolve(eq3)
solution_eq









    Out[32]:





$\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)}$



In [33]:

    
# Solution

general = solution_eq.rhs
general









    Out[33]:





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



In [34]:

    
# Solution

at_0 = general.subs(t, 0)



In [35]:

    
# Solution

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









    Out[35]:





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



In [36]:

    
# Solution

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









    Out[36]:





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

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

df(t) / dt = r f(t) (1 - f(t)/K)

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



In [ ]: